<!DOCTYPE html>

<html>

<head>

<meta charset="utf-8" />
<meta name="generator" content="pandoc" />
<meta http-equiv="X-UA-Compatible" content="IE=EDGE" />

<meta name="viewport" content="width=device-width, initial-scale=1" />

<meta name="author" content="庄亮亮" />


<title>Bayesian Linear Regression</title>



<style type="text/css">code{white-space: pre;}</style>
<style type="text/css" data-origin="pandoc">
a.sourceLine { display: inline-block; line-height: 1.25; }
a.sourceLine { pointer-events: none; color: inherit; text-decoration: inherit; }
a.sourceLine:empty { height: 1.2em; }
.sourceCode { overflow: visible; }
code.sourceCode { white-space: pre; position: relative; }
div.sourceCode { margin: 1em 0; }
pre.sourceCode { margin: 0; }
@media screen {
div.sourceCode { overflow: auto; }
}
@media print {
code.sourceCode { white-space: pre-wrap; }
a.sourceLine { text-indent: -1em; padding-left: 1em; }
}
pre.numberSource a.sourceLine
  { position: relative; left: -4em; }
pre.numberSource a.sourceLine::before
  { content: attr(title);
    position: relative; left: -1em; text-align: right; vertical-align: baseline;
    border: none; pointer-events: all; display: inline-block;
    -webkit-touch-callout: none; -webkit-user-select: none;
    -khtml-user-select: none; -moz-user-select: none;
    -ms-user-select: none; user-select: none;
    padding: 0 4px; width: 4em;
    color: #aaaaaa;
  }
pre.numberSource { margin-left: 3em; border-left: 1px solid #aaaaaa;  padding-left: 4px; }
div.sourceCode
  {  }
@media screen {
a.sourceLine::before { text-decoration: underline; }
}
code span.al { color: #ff0000; font-weight: bold; } /* Alert */
code span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code span.at { color: #7d9029; } /* Attribute */
code span.bn { color: #40a070; } /* BaseN */
code span.bu { } /* BuiltIn */
code span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code span.ch { color: #4070a0; } /* Char */
code span.cn { color: #880000; } /* Constant */
code span.co { color: #60a0b0; font-style: italic; } /* Comment */
code span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code span.do { color: #ba2121; font-style: italic; } /* Documentation */
code span.dt { color: #902000; } /* DataType */
code span.dv { color: #40a070; } /* DecVal */
code span.er { color: #ff0000; font-weight: bold; } /* Error */
code span.ex { } /* Extension */
code span.fl { color: #40a070; } /* Float */
code span.fu { color: #06287e; } /* Function */
code span.im { } /* Import */
code span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
code span.kw { color: #007020; font-weight: bold; } /* Keyword */
code span.op { color: #666666; } /* Operator */
code span.ot { color: #007020; } /* Other */
code span.pp { color: #bc7a00; } /* Preprocessor */
code span.sc { color: #4070a0; } /* SpecialChar */
code span.ss { color: #bb6688; } /* SpecialString */
code span.st { color: #4070a0; } /* String */
code span.va { color: #19177c; } /* Variable */
code span.vs { color: #4070a0; } /* VerbatimString */
code span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */

/* A workaround for https://github.com/jgm/pandoc/issues/4278 */
a.sourceLine {
  pointer-events: auto;
}

</style>
<script>
// apply pandoc div.sourceCode style to pre.sourceCode instead
(function() {
  var sheets = document.styleSheets;
  for (var i = 0; i < sheets.length; i++) {
    if (sheets[i].ownerNode.dataset["origin"] !== "pandoc") continue;
    try { var rules = sheets[i].cssRules; } catch (e) { continue; }
    for (var j = 0; j < rules.length; j++) {
      var rule = rules[j];
      // check if there is a div.sourceCode rule
      if (rule.type !== rule.STYLE_RULE || rule.selectorText !== "div.sourceCode") continue;
      var style = rule.style.cssText;
      // check if color or background-color is set
      if (rule.style.color === '' && rule.style.backgroundColor === '') continue;
      // replace div.sourceCode by a pre.sourceCode rule
      sheets[i].deleteRule(j);
      sheets[i].insertRule('pre.sourceCode{' + style + '}', j);
    }
  }
})();
</script>



<style type="text/css">@font-face{font-family:"Open Sans";font-style:normal;font-weight:400;src:local("Open Sans"),local("OpenSans"),url(data:application/font-woff;base64,d09GRgABAAAAAE8YABIAAAAAhWwAAQABAAAAAAAAAAAAAAAAAAAAAAAAAABHREVGAAABlAAAABYAAAAWABAA3UdQT1MAAAGsAAAADAAAAAwAFQAKR1NVQgAAAbgAAABZAAAAdN3O3ptPUy8yAAACFAAAAF8AAABgoT6eyWNtYXAAAAJ0AAAAmAAAAMyvDbOdY3Z0IAAAAwwAAABZAAAAog9NGKRmcGdtAAADaAAABJsAAAe0fmG2EWdhc3AAAAgEAAAAEAAAABAAFQAjZ2x5ZgAACBQAADWFAABReBn1yj5oZWFkAAA9nAAAADYAAAA293bipmhoZWEAAD3UAAAAHwAAACQNzAapaG10eAAAPfQAAAIIAAADbLTLWYhrZXJuAAA//AAAChcAAB6Qo+uk42xvY2EAAEoUAAABuQAAAbz3ewp/bWF4cAAAS9AAAAAgAAAAIAJ2AgpuYW1lAABL8AAAAKwAAAEyFNwvSnBvc3QAAEycAAABhgAAAiiYDmoRcHJlcAAATiQAAADyAAABCUO3lqQAAQAAAAwAAAAAAAAAAgABAAAA3AABAAAAAQAAAAoACgAKAAB4AR3HNcJBAQDA8d+rLzDatEXOrqDd4S2ayUX1beTyDwEyyrqCbXrY+xPD8ylAsF0tUn/4nlj89Z9A7+tETl5RXdNNZGDm+vXYXWjgLDRzEhoLBAYv0/0NHAAAAHgBY2Bm2cY4gYGVgYN1FqsxAwOjPIRmvsiQxviRg4mJm42NmZWFiYnlAQPTewcGhWgGBgYNBiAwdAx2ZgAK/P/LJv9PhKGFo5cpQoGBcT5IjsWDdRuQUmBgBgD40BA5AHgBY2BgYGRgBmIGBh4GFoYDQFqHQYGBBcjzYPBkqGM4zXCe4T+jIWMw0zGmW0x3FEQUpBTkFJQU1BSsFFwUShTWKAn9/w/UpQBU7cWwgOEMwwWg6iCoamEFCQUZsGpLhOr/jxn6/z/6f5CB9//e/z3/c/7++vv877MHGx6sfbDmwcoHyx5MedD9IOGByr39QHeRAABARzfieAFjE2EQZ/Bj3QYkS1m3sZ5lQAEsHgwiDBMZGP6/AfEQ5D8REAnUJfxnyv+3/1r/v/q3Eigi8W8PA1mAA0J1MzQy3GWYwdDP0Mcwk6GDoZGRn6ELAE09H/8AAAB4AXVUR3fbxhPfhRqr/6Cr3h8pi4wpN9K9V4QEYCrq7b2F0gC1R+XkS3rjKWXlfJeBfaF88jH1M6TfoqNzdWaXxZ0NM7/ftJ2ZpXfzzeVILi0uzM/NzkxPTU68Md64GQZ+vfa6d+P6tatXLl+6eOH8uVMnTxyvVg4fGisfhNfcV0f3luz/7Srmc9nMyPDQ4IDFWUUgjwMcKItSmEAASaNaEcFo069WAghjFIlAegyOQaNhIEhQxALHEqIeg2P0yHLjKUuvY+n1LbktrrKrOgUI/MUH0ebLc5Lk73yIBO4YeUrL5GGUIimuSx6mKl2tCDD8oKmCmGrkaT5Xh/p6rlphaS5PYp4kPAy3Un74OjeCdTi4nFosU6Qg+qRBsoazczLwHdeNqpVx3AW+oVjdhMThOo6YkGJTl862RFq5r263bbYSHyuswVrylsSBhHzVQKDU11g6hkfAxyOf/DVKJ1/HCvgBHtNRJ+b7eSYepeQ4VLZBqAeMjgM7/zyJJF1kuGw/YFpEq458Xrr65YTUa6VCEKGKVdJ+2FoBYYNKCwV1K6B2s1mJnPB7Ww6GtyO04ya/HHWPHs5P4J65NyVa5VA0E0LocwPci45b6tvMvohm1BYc1h12Xd2GrbbHVkjB1pzs6IKtOHeYd+JYhFasmfs9Zt+SZlo9pu8eg0utWZAKB8vjaxBQx7cSbK3Qdr2nBwM27vrXcUHtLolLJyJjK3CAbDcFDo3hsPZ63IH2RrsoWyskdB47jiKitFtcAgqj4wQQxN3PB81RCiCo0Y1jnUVYlOj5JHhJd2JBevIEeSQxDWzTN8PEE3AL90KtP11dVrC5II1L1w331pHFq10vPBGYeyUCFRvB7PAEzMltdubhb+lZ4dw9w86yyNfG++u0ZWOBkmsb+GrsrKGIN4R0XPQimnAEcj3CI6ZDR35zzHJEZlcW5cQCTMwty4umkB5B4ajHwVNhQDqdMLSAmClnhLScgYgMbQJESALUrtIvjpQz9LVxuIPSiYgQkjusZ01l4BERrPtdO9KfDErKQLne6EUbJlXHqTccNzL163tuES26ickjo5va6FIkCyIyaFEYA+lejuqlFxLWIYKmQG9W0tlMe0yXu80wPe/OavEJrd8srSFziSal30wMj5H2mH7T6H218RQ93qOFysDEgtLBoRuQUeXjyPQKexdLjoa4vtAQJiBsEXYutEo9T1/m5mUdBMbXFCzIq8Z6Yl5+7nyic+1mE3xisVatpBarpcC/mUs9/s3Csty2GRPfLMo7FrfqcS1KDxIntwVjnkEtjRJoFKEVHWmelIyxd7Y9xlqGHTSA0VfbnBks08M4W21bHczuJBrTiYixiBnsMF7PepCwTAdrGcy8UqZb5uWGvIyX9QpW0XJSrqE7hNzjjGU5u1vgRe6k5DVv4DZvpVnP6Vi0yMKLOhUvPUq9tCzvFhi5mV9KVNMvWpfRJg1bggjEml6Uz6KmiiN92dh+Gg19OHK4TmOC61TIcAFzsF7DPNQ0fkPjNzr4sMZHaEX5fk7uLZr9LHK9AW9KF2wU///BUfaOnlREfyrK/rv6Hyn3ISkAAAEAAwAIAAoADQAH//8AD3gBhXwHfFRV1vg5974yvZdMQspkSIYkQkgmhdAyIIQQWsSADCLSpajUiMgiAkuJNGmhKyJGDCyybCiyiGBHRGQtyLIuf2UX19UPy7oWyFz+972ZBxOE72N+L2+Yd+be0+5p99wBAscBBIN4ACjI4D4oUJEIVAbIL8wPYX4oP1TQ3um3+0v5dZz2bj44nsyKLhYPXKkaL1wCAhuuXcQ69dsWyAu7qF5PBMFqQzQRkzQgYvIQCuXleXYHlCXl2x1YZg+F7HxMDNAQLQoVetwuKZCZjRUTQqc/f7RjebisqAeuEQJXmpZUdA/3KgcgsJA2kL1xDNPDZqCyQAWdXiIy5YOHThUq4/KB1XFpgPr5heVtJuSQvJzxOeKB6HfEplzKWCEA4Sc+Vgqkw8bwIF16K7fg0ttNJr3DajEKBqfT5UlNkwXJKyD4hCRRlFySwU+TvTTJkJTh1wkms6l/pBWa08Fmt/WP+Nz2AWYcYEez3WwXvU5qECE/VB5ylJXl5993Hyc3zw6hkHaPoerldxVjh7eMX/F3hYWxu0KF382pcKpXsV+9QlS93Mj/Sz/ujinsVE1dDTszcEk1u4LpPdjXmDdw6UAsqFlUg7rmf2J+d3aGLmC757GBuEe55mHNXGxifZVrLtuNNUBhwbU6wSQ5IAOyoS2MCxcH7VmpXkHIdZlFP4BPtOvFdvlZZsncL0Kl1pZcS99Iam5eK1erfhFvrkviL9HDKc5X6OV/ChUq7aGEvw5U6QuFVCbEhOSSZHegODM7WOzxhOzZ2cVFJaXFIbfHK2cH7WlELuK3EnR5vHZJEkzvHZw35S933n0ucur5ky/MO7SraN2mrVuqGiNPnIt+NnTy6HF4fMkfvf+6EEjfkpWPh7rtXrJgp+NAk9hzQScj6194/+yxlZE72Ow0KvcdloMLbPcBiDD+2jdSW/Ek6MENfk55AfQMtwabaPC0aZWZ2a6Nob1NKgxRc3qemb/aF0jtk3xZPtkpc4Xjr3KVXE7WDfpi+sfVJ1RotwUyJVFVbE4ZV3JUPi0pLsq++XMM4A9Vd+/YcXcVvrtx7bLN61av2oINVTU11dU1NVV4cuPaFRvXrV7xDGPNH6+heQJpbMQaHLiz8R9fXb5w8dLl5vO7XnzhD7uef37Xxa8u//3ipa9pxpUqrt5AYeq1b8QPxVNg5BQWw13h9k4PpEqB3Lx2eW0DlmxfqkdfUhoy9Y6EnNZgW0t7MZ/6smlubka+I0NfFckQoDwPkjih+d4yrpTleTdRqoinJE6Ts7AULcTt8mRxQbYjMeLcXMpYwucgMgaCkrrMn668Z97YBwZHJm/+/hnWZ/KwOzazl5c2DerS+o2Xth9eshXXd7jTu7NHHeb98+VHfqw/+z/Cmp5zhvSZe3e/kSOubt2EO3tExnWrrbsy/51x94+aWFa/84V1k/bfx2Z1fWE0+2It+2zfxGEfAaBiMbBctRiug0CpIBLFUpyK2R+OumYgYrZB+cZAdoT4+TfM0CpsksEggGCxGoNUsV4J5sVpc5SGJE6pwxvIJgM3r97+1Kq1S7et2UQKUI/v7znOCn/8jpW80ohvKaN24aOatFEFAx8XLFYDFYItR0UbkQMljuIiEgx5HMS0efW2pWtXPbVdGZb9yjruPIInv/sR3z/+EisAhMFkrmCRXGCB9uEUKgoomw16o95qEwxoJiaT2cDtl84CUP5G4XWJOTBmWLK8olOmNOjMKhUpWZWHK5LZgl9279229we2OBUX50kuVjv5QDo7PBwnsvrhWJF+YDIuVagZDxeFHOF1MEKbsBMEQS+KJjOVdXJ1BKw61EH+feqSTzTz3I7ZA3Zuv+whshy3sDFL2TjctJR6n2SDsfFJ3A0I5ewXfAgugw7s+0XQG0SAfFVWHOEsr6TyphSHW5NHFc9J6Wa+7B3Dfp42HguHAUINniPlZCpQ/l0CogDIrW/8u85iv7sGv8ZzGzYAxjwV/MCxTwobJQCTWU8HRPQeruaaXpRqestVdUOXso7dupeF7px4Z8+ed3arKFc44AIg51W9ch4kIIiUEocmSk4sBpCcj15oUDRJXYYExl37RmirrkIv55rLASYJJF+S3t0nopeptU+E+mLrLK+lPgQyid3mCBU6UP1rVz8R2n770zc/Xf7x8s/Nn9fvaFi3rmFHPfmMLWRP4lycho/jNPY4W82Os88wiJ34K4tdAIQjAOQkx8YArcM2PaAOjSZBL8uolzAJFFvGDXd8ej67P2AvKpUkOYghcnK7zl300RBcsExwzJ/hbrd7GuYBwhgAIYtbTx/3+d4klJ3gtKCQnGIz9InYZEzqG8EkjSzNavCB/cXYlcQshhyMsZrI6PYLWc3lOG/vlA4rHr/3uTFD3r38/r+3fMKOke9W4oJ9G566u7au84CpOz/ct5R99wF7W6dIYjjnawrHIAh3hlungFOWgXoyzVKbHOr1eD19Il6vISsrrU8kSzbY+0QMGpdjgYh60zDTHJKHoyP4404pw27zB4o1o62gq+BLL299am8j+zv774zj995/dgTOZsOfWr3rnTWPj2h8qGbo1/M//kYYvmxfms7TtPrM54E7ns4vwBw0rFy/aNJjRRVTet31OgCBPABhongUDOCAzuE0h6gnxChToCJ1ulB0iH0jeqvscFBZotflk+hMQ5oJDqhrC/l//FxmAUlGYeK5Z6Jl5MDec2yJQdc+l5ViNduL1avoZ805eGll04jy6COKheT8S+U6kQwdw+lW6nPpXF4qtEoBziwAye3mMnRLkqlPRLqZdQlsKxTcLghkqhzjrLL5M+WgUwldSkjbL1HPLrCf51d8MHbv66zu/mcGl5Kz0YNZ0+mcf759kbEB29qGGrZiYWop2b2R9fYqnKnlWOVzqXqgNfQIB5LtRr8fQLLT7CyT0ZLaL2K0WFzU5e0TcfmojkckcgvcyhJ4pNlr8Bd63VyEhIbiGhfIBFGTq8R9lqcWB2Dl1G79Rn/9i8n08OU3L/760UX2E369YuvqVUPrI9VryFR8CXc5V/rYefbW7svv/YNdxUHv/OnFVQ1V8yse2Dde0UcAIY/zU4L0sA1FEQg3jJT0jVAJFBlqbOOrALk1dCOmkuHNF+mpaKOYunHhldNAlZhEyFGpz4R20C+c47Vmu+6gqXo9lewuq5TfXrLnZORk9Ink5JjAlNwvYvJBoF8E5N8qd9nN3jrmj7mOx8OPLDXqolpgwv0zZkpuzaeTynf+vWjNvnr22b+bsfDJR7+e+cL6dQ1bXlu3CDvOWfHIMytnrhJPHt7x4L7eg/48+8C5U0euLuu/f8ozr1xteHTRssdGru8V3kwfeHTMsN937/zksLEzFdlO5NQpNsMLWdAtnJlizzQYAAQu26AljUvWZbEQlyuJi1Ymcr8Iaal2jjKNg5qJ9Ctqx02jMyDFKHJw8TpUIvjHKhXZQlZ0/Iwe1eO++6/RVHpg2mv/uPbBuguPMtfKLU+tuXfjkIFraEVzg2tlMuZg6O57/vXBP1C3kZ3H9od2PPV81RMVE/aNAy3HEcaokRS34Ta+LAA8XotzQMRiizkRDVfN87X0JXae6NzkVR6Znehb6J8XL+Y3IKovXMjn0oEDMrkmmc2iXu9yGm0DIkab6hgTZklwj/T6FDccpXsmn6Rjlxv+knyrTFMR8+U/cF9+DiRwh/UCiChwdeXD58cDhSwsRjeikNNcTo83/0AtP2DDKLywji1nhxSezMTjgo9eVHOy3LBbJgIQ0OsEsToiIFRHrIjI4wHOlfxEz6a4ZOTXTLq9eTjdTofW1bEH6up+g5GIBDhGEr2BkRNVlMZTa/P3HKVyrMMKrF3H/KPYUAWjlGsXaRnXrxTIhrJwqp/bMtnphFYWIdgGoLWtddqASGuPzdA7YhNaqFZLvVJSEa48LZwUd4YSN4mJ+aq/ctSSXgtmD6gf2emV91/9KNj38bHd9l3PX0tq19dMnzFw3OSsgsWjj+zqPXn0w4On3e9nZ+NJLYFZ1yqkQ2ITFEM5zzwyA+1KLJ1kVwpAjsvSTgx3S+rQQeiisxv5Ky+9kGbnqUmllmSFEhOP6/G4ug6C2nJQUPdSt0td36R1IFMgbsUalrqlQAbw4KK1v1BwIH/udKqm8NCQbeMHP2LUtVk3rv7Fb4712N3Tt/DeaWvZt3+8wA7swe6Y/5cvjv3I1rHJn+AyhLM44ODVn14/7bBUDpq/hpxb8c388XfdM+rU3veu+Tws17Pv7O79aFvzMnvxc3aaHRq8sAZX4jgUsP7CfvYntoNhGYquJiAAAKJNPAIyWLjk0ojFqENR0SwqyILNaiG9I0bRYhFECoKD518xh6iplZYz+5W8H0OIlBsz/tURB6IHmnaT7itJORvb6A94cnbjGZYvHrnSg0zENwfPGTGddQIKJwCEo9xyW8ALGdA7nO0UUg1Wn89iEGQLjwd01iRrUlXEarWAxVcVsTjAWxUBevt4QnM9/gxBMbluwe4SAjxpj/mcgN0ef3cCt2IAhVVLsR/7+TIjjZjU9PTeY1ew4I9/Ovhn8cCeI/Nf9BnK2Pk3/kZ7TF00+6HoquhndauXPAGAMIdb09Oqr8gOu6jFpbdQb5IDekccglHi/HK2DL+4emRymUNIE3+Ro3WokKfbtNP37Cs0/7rxjQ0X2Cvs2Rex/NNLuysbxBB7lX3FPmdvl64rwyU44QusOVSzuj8AUTgmDuEc04FdsYcWQQ8COJyiuSoiUsFSFREct4ppwc9rSBlA+ZuAPZTBx2Az2Uo2CY/hIHysic/1z59PI/dU5CtWz+aJB9gi9gKmYebVKZgHgMq89Bc+r1GJWSSDAQXQoWAyS/reEUlCQsTeEUKRr3B03DZmUZBwxy/6S/MZmh+dTYZHt5OF4oH1LKc+eilhJj0UhpMlAKQ6pAbjTRPxSW45Q0CbAac3asPzwaNfrY9LTuyi2ilOhUvnI8SSohNapUJK7wiAaDLZe0dMgujtHRGdt4+8/HaphRyV9+rq5lT1xe9nfPc0a2IrDuKQL//9bve3DrL/so/Qj0kbVrGXCYuWZWXjUhzzD7xn/+D6GvYau8Q+Ze8H8LUY7WK6yuVQ2KdHBJ0giCCaTTraO6LTiQaJoshJV81RgnG/Qbydi5f/DYnpjc2ssZGSRrI3Ws1z7dXkYQC8NoLNxfFqVpwaNht1OotVT4GzFDJj9GrpGI15+JJiPpxLMg0v6dVv9AONx9jclFWuR6fyFGvI0TNxvRC+UjHmnkjBViRGg4Ix0Yn6RGzLWkgJZRVRDKHw1TvRrzc2NpL1J6JN5M0l0dc5snnk4+jCBF0QIT1soQCCJCMFzgtw3EBXxTekkO0+0aio0pV/bIp9V+KIgpPrUZJOFCUev/JSmsuNBjuVjDK1gKQgp2DnLbuZlRjwuJUAn2MY4nce4COtZjadZSsCntbhh6zRomMm0bbpo+bh4oGrVQLPOume7Uev/BCXo1IDsUG7sFsvcaytVpDB7jBS2aqjKCdypaUI4xPzabNJKZdj+WvNn+tsW4/RVB2xkGeEk582NR/nE3ZMwaxy2guAqFp99FZ5bu+IXqDW3hHqvLVNiOltBiTmueJRtpW9oZgjHIE9sBOOujo9+v1/fvn5h/9Eeb77LHuYa+94HIt1bArbxs6yU1iIuRjEAnYqZp+E8erqdUBRONnA+c75DE6XQaiKGAySLDuqIjKVEtavhpXmSgW/mlplYChutYXx7Ay7tLsRZ5PWUePGL949euKoYPr7t1HOh2jK6mdXrVC5wHaoXLBCCp+Zp8MeAIEa+OqmZtns6x0xC7KTL2yZM+MtlRs3J6I2pViG8q258sX7OOxndrH0tpz5ki3rzuqxivyf/DnN+WMCN1SGs8yIxKS3y0aDQdYTwePVm8EMVRGzmVDK5UepkSi6cntnp2Ku8ktw20SOf5bGNm4BcRXyGdhfcfkJ9jQ7/VXTzl2vfEZGRLeJB94/zf4+LjqZjFi9cuWqJwDVHIFw29ha4V6a0wSQ5BSFrGxTGvV4uH30CFSfoEoJiY4mt0CGlozy8D+o5jgx+6jmBbwy4BEI+9d3rHnZ0I/GN+7usnL1ey+xM389WLx/1+INHRbWXfoDLjz+6Z07su+YN73vyIFFvd959sV3qtf2nfFA35F3FQw8AoDgABCGcv7JvJ7iABSRUp1epgK3CYLmFeJ5qGYSi7k3IEsbWYFQyQrE9PWqJzjM14yPj2OHrLDdhgYZZafDrqOCmQ8UpzGUuFzsLkUnVHMYs4uij/2F/cJfFxrfee3ld8QDzf2vsC8wo5nuaa44+Mabh+ghQAAA4XW1/pMcNqJgMuooCJQqiPLlrxWvQhjgF8//SgXTwej3O6M/NmF1x8zWHdVaFh/5uU3bnwXkmg1yXz6aT6km+QwpyW6LRdQn2Q0U9TGTotqUGOKqNclWAjJldKcyenwSZ0h8cyc75y5CT3v2xU42u+nL9p6UYpSa0Nne7yy+1EQ/7PaW6/dbm0N88llHNx18ic5qnrv59RXv0YUK93QAQr1q9QNhhyCJ3ORLiskXFJMvtDT5KhocAz63Yu7rj/PIY0oTXmKdjuAkfHg/60QWROeQZnI4+gq5M9oX4lybrUY5GWGrIBJRpnoDiChTUeOcJmE+qKL+GCJdcNEhlrSb+Q6T8+R887zoCZJPFyv1ZQBBscZ6pWKmQyqDLKBgMIoCNwcUdUrMcuuKmVot8AvlzU6qi9roq82/0LSFwoaNC69OAIQGdoRMVnSRY2mRUFAYoxcJlTDIOdBSfeJRD5nMSvEEu4B+dkS6svyKX6HWC0A+i1c2Kd5c2XRy3h0mgYbo/4spg/KNEDuCzdrMFFACSacHOUgFevPMXj5rMb9CfMoLfOrSA+KF5b9KyigFJCgExOMgQVJYD1TWiQQEwrO+G5rpVFUTC3DfaPxsA1vG9pEg3dQ8jnwV9QJea2Zv0k3XKtUKsJLHIlEqwBgjmU/LQUfRp9mbCwCxTjhHHZIf9OA8AILRID2BkJ+s1ZoxwDW1OMStBHU83G1fm5MZ0+4QzhUdK3f33F8MRKk50lPCUEXzoVc4K1NnTEvz+Rw6yqMpYkzrFSFGI7jd1ooIt4LJFRHRA24o/98LVH4tX7NllapJZ7zS6LZn8QVeLKsVKjrQrxv43GPPvUychyc/VveH0F3HR77xCrNs/mPDWy89tOWB3js3Y1+b1GPe7Jq5dxTuORZ11TZuHC3LD00fOhwI7OVWtVZygRPSeVUt0+D1Wq2mVGqiGX4zmNwOu8HOhccRljzgqoiArYV5DSXF1SDB1sddEk825YBijeRQiVcrvHAqyJ5Pv/3+k0l/7GwKzGzQ6Wa811i/qXFjfb0wlJ1jP/DXxwMGLpdcbNHcsTuWvv7ll29fOPPJXwAQpnMOLxWGxbIaK6VuPU3ySmaOmQ0cHDPPzVmNGM9qlJ1DHgNzu6hmOGTcZXYV9f8d8HTbUOn8QrbvuW11Tz3swiw0oRPvyPQu96Sywe9+2mlNGRBlVqGU88fB+dM97E+VvGCx2CV7ht/htgIgmqhez9mjt1FnRYR6bscerSYTkLTqvTcUDPLPA6osi+JOiG7ST//n2W+/++TCTLMsNCxmTzdu3Ny4evOmNS9gNlr5647tA/rh0V+/mfny+4Gv3r54+i+fxLF0cN44IRk6hdOTDF4jpdzqtkrxGit4uRskyaUyyqIw6paZQyiRZQ632++JsUuivNbh53Kb+x/2JYp/e/+7qFl8eecf/zBk65bfb7WQLstc2AZl1GMH9v3fJxx/p2pttp/+c/eGrS8oUksFoBYpHVxK3cVlMjkJ4UaSuj0GvhQMgKIsVkScspUqq0GtY98IAxWmOZS1p2QNgeJSXkPW3DX3mE+zrxreeANH3lObN6LH8KHopW83l9G3+3TugmsDC9PnPNkLgEKQuYQCzplcKIVu8HC4a56vQ5YpvYtY4ESnSHIzW6Vn+Qzd72xlLbYWV0R0nXpFDJm6XKvOqvPk5pJekVxrm/JekTY2T7teEU9KnHUa+zj/8pXd+rzbxD1uragaVBdAqDC+jaAUkrJv/OXKcGMXmJOnbhQXF/F3QsHJVnf87VhB3sSqoa/te5X9jf3r7FdPzMgtC/ccNOnTtwb3ZPb6ZWdOPLzh7amPD50/4z8/1T4uVE5ICkzt9ewxXYdBbfPqVx54ddvqMauTndXFnYfmBnY+2PS66ypEhs2ZFOn5IO08/ZFvfn4cEPYCCD24nnuUzM5i0nFz7dF7vEkWvcMhVEQcNgOA3q0Y7xjlCatesVT2mALbtRUfM1P06cfm/+GZhgadoWD/jBMnyJuLfn/kk+jrfHXnDOow4N5XP4gWAxDYDoDjxAtAwcr9tZ3PJCDa7Ga5MmImVlQ04/3EwqZSIqAJJVQc3NDQ1CG3TceObXI7CJWYU1Zc0qFDaSkAubaKudSxTZAEd4Q9TqPRrNP5kj22yognrLcC1z6ISzW5xSTOhATTljhb3v2det7Zv/eNGZnLt9g16B6h+aqNHZHv0yaP8TSV89QGJTzetxgMRqNOEkSdYHeYAGw2nY7KRje1xiKGfD5zeUyFyuJsRTUiQi0bdclYkzcER73JeuD5E2zOnB07dKSgy2icydpGlxLpQTZOcjW/XTo9NjcO5nNT4GQCoiASQHfca2tMVBjHYVRo6SRfJQGoCAfcdruDiz+gdwRo66xWHrfb4RPMPm5p0302p1UPDkUPuCLEt534Igi1bHVIVIgEzfAqepHh1bRDypryyOa1DVNmblnVsDhFl79rIuIAXcHhmYdfJicWLNj3cnSLcv/zx9HjQmV99dDDg8e8+heuMZq2cnxdUBBOApeiri69x23S22xcWW02g/V2ytpSV72Jmrp7m4JG6NDUt95RNPXwJ+q8d0XUSWM2dhSfU9EknsU6wSyDnOwzeLgds1GbYvxvmcVylSHFilGFxE4PYRT74fKaf/wOTZcvobX5lZ3PPffii88/10Cy2I/swyeR/AFNmMfeZ1f/8rfzH545p1j5vdyW1apU+6E8nOEzCrKsS3foHJkBwQhWq7siYrXprboUaHXDzMdZ0GLBqpaeO2hPAhMUr62Y+gRHrThpU8Niry7c+PBf/+f7yzvryabGFc8+6xowcMRg1kUqqh9azT5h/1GcNr14+GTWl29fevfUeYVXHNNSlVexqMKW6qHJyT6bL8OfnOK1pqalecxOp8wtv80MFRHz/+Y2VT5yJ1l63Ul6r3vQ0njtQyL9GzaIW15cvXnjnI8uf/fJ57P0SQsajObpM/d9mHXp3YunT59birloRDO2a6z/9T38eEzFCzE9okGOpw1ywy6zXm8wEF4DsZrB4FYtg03rc2nRkaE5IY15ZEfvjt4eRQtfaahz6rrsFoaZNlk/fTbaJFSenDQjlrnS6XyW1twOtIplrqLzeuZaEfHYJKq/rj/5t8pdueG5kbsG25Hfpq50+j/e/+tjA/bXzF82+dmN88r/evSPL3Z6ftEjj7Yds+J13jSzsaHnpjbt7h4Uvrdr2aAH+yzaXLm4R1W3O7p2KO71FCCkX/uG7BQrwKPWJlwu3jPioEKS1+C0OXtFLGGbVeaCkj1xU3kqIVjV5ONWqo52xVGXhtxKNuHyEMcdA5NSJuSy17ZurRiBXdlrw2vN8lyzHQeQZdU9/83mRWePngiAsIOvrjKhElx8fh86ZZPJ4DS4PSaz2aZzWdVV7TFqEbMS/4daVmW0rJcrhBY127EvX9TPNNQl6UP7Z7zztlAZLeMO6GMSvnpozV2Dj54hp7RcjgiVau+HAQ0ms6hHK6jhiJZl+NX0NFTicIYQt7ER+76ptuiMte/tYyP4oI/8o0cx9iPtrx6K5UpSgI/Winsblz4lNc3rsZipYBZ0yQ7ubnTuxCyYK7c2A1U2Z2Rlk8LhUHSq1BmbsoRPKeSfcBbp2qSdPsY+3jNxsk5nLHCcaHqjg0snBF7dzc6QBZ3OvHR/dK5QyUaz6j5l+4tJbXTp7trW9eRvHClACAIIOpXGzLBdFiVAUWlxQZ3RLaD1pnQ4ngmjmhUfYgteQT9m/JktwFVH2Cn27hFSQLxsGO6IfhU9jUdYD0AgfL1LfHw3z/sVMqnHK5jB7OBLO0UHfIJCVam1GRJo46KKOdrSUrLvuwFOnfnuS/tYTsWfl/StKu2xq3cXzuCVn9wf+pn87mrGy5vtC03HtkAsZ6YPCZW3yJl7RUQr6npF0P2/5cz0oeZ/ksHR0+TL6D5y31Q6eN685sPxrixetlPl5/YlJxu9AFbZRbmnpqlpTq09K3F7TdV/bpXcPJZTfEtxCddDvj7d3EK4ZLfHjedrpx794PFH58/49MClCxdM44aRZaRxE+aPjywnw0Zg4ebdS6Xj7NzZoCl4FhAvMxuZrfluorSo0RSABN+tlHzx8nKeJv3cDAiV7Ijaw5Oq4OwWDQ4H8UFqqsXiE2laujso0QScEzYFFXSDxYr7U7DPVNCV5Dj2pcRw4eKhDx+Z/9jjp45OnvHwVFIePIvB49LSPRvZ+yPvJcsjvOq5cRenZNg4zJn2qEvdpyXVQg6tAS/XAzu1JvkcpuoIdVglCaojEuTngS3pjfw38rSkOlOZT8nQVNOmbD9lKoU5HFg8t2TMUz2mRrqPyi95omTcisrHK/sMJSfuLFn/UKvsVinhsvqH/RkZSeoOPFuKdcJwrcuYCALV8343AGpSu4xtNPOWXcZcCQNO1/Xt0PNKk/Gszp3Ly0IVZPfVC2Lfxb3C5ZVhQDjK7fd5dVemazjNozNTahCARxo62irVJxKnwUz4SzDKgg+07k9ljt9sw2apra1KOJCldLR6NAOuqD89OWHNwpPHcdniPisKChY+tHv7My8sX/FdifTO+xlov4LNXXfvoH7vstCH5z462QkQypUYSDzBpV4Zzk5y6s3mZI+dGD1OMS3dlORL6h/R+3xOcNr6RpxJIPa5uRWkRdPQzZ6Nm29lf5Lfinl2ypuduEqQxqONXTatnD0HG9jQblU05erVU2+99f/EEzUL+/1uGTs397MxS+7YtDz/xwtzsfO+U4psZqMkeIVtnHNByAibW0GmBSxtctLd7iwZeNSYn1gJchaVBku9il8r9co82Ja9clCxDnKwNLs0IXQ6VLV4+OLx8+eOq7t/UVXVgmF14+YuGrN42MKqeVtnzHh627QZW8mHj01aNmxh794Lhz059ZEFD/CHvfj7JZN+N2XbM1Onbd8BiscDEJT9Fw8MDrdzWGSj0WYS9URPTS6LW/YmGSwW2So5HBScbqsz3UmsTqvThG7JlATlWg+33RHrzL7lpjuGUOGj1uaovjBEKnH2HjYCJfY6dmGv72BvYGd+ARu7j1wgZ5vZ3Ma57Ec08RslQBKsgaxUVYkkUR726QUqUDlmFjgmiYqtbgjFLYRiI5p/YebmnxVpXPuF1kupUABdeGdcdiE4pdy0Dj5fmkmCgNS13E07lbRqK/n1/mCviN+tt/WK6OGGznh/s4t9I39VVFmLztSUlwuwZdCiRC2l/Kk33lG0dHD/qprTbw5/ZmTxqMV9Z8yYvelw/cCqjf/+6K9P9H9t4KLl7R+cvmJR99W/f6Ggbs3LPQbRnMF1WW0mD5q1NDW4IJjSKdy5prTH+klDl+fctXrZxm5rs9r27dWuY8e8oqHTRvWb0MVZPfnuKWXOMUCwWLTQ8eKH6u5TWpiTanKAI8lnpW495N90QCAhzctKeI/FxVnZpaXZWcU4pzgrq7Q0K6tYnFrUrl1RYUFBYfwOQGEM7xzvEdt5hxKeSwWDXmrNT0936a1esbSDZAKH1ZRuIuCwOYjJYXKk5AWcoRQByhNPBdhblgFRMxHuG90bnN2obu8KDjc3eYHM1py5DiFU2NqhNXTQOXMWz10weE77sRWvffDZq0880vHB5vXv4PB3les1tv2D02z76xP2YNvdezD3pT3s7N497JOXhMCeTTu3t/2dq9X3n575qfMjIXZI/Q7b/u6brOGD0zj0rT+wD/+wB3P2xr8GQKCCushU8W1OdzqUhlt5pRQDokeJazP8rQwGh88D1EYJNTvSOakf3feGku9qVGpqG4xTV8ojfbXWGSt18iYUtdZJXEnDlt0/edPztWvHjM+btnB+HauecmLUlAeov2bk6HHjJkhCcGFoRIcJs1jnI2OaCgRBqd8NhFraSI+CBGbICTupxI21YNTrBbMkWKwmUYegHGS5WbPRiyhjVuw2EAfPVEriM1kjLsUhtexzTK9lO0kQ1/dk29mzvXB9yo23qh9EHfeDXhAhJWwiKKAki0J1RCSQr20nattixUJOXfM71Bv9Hhc+CdeuaV3LRAIbAAjXdUoX16r7wqGgF3iOLui5Zpn1JodXKu1gsnFoi9Pi0DmtjnQHAR63E4fT4bythikCCP22ZKVVoUS+hp0Bqm51Fnr+L2UjHz5YPXLwfRNx36B+l3eeXrwWxYbNVy/8n+pGrtwd7tNtSfXsNFaLo9jTdPZ89ub/pXB47YrkEiRpzW3r+oJ09UfBJLnmAoG5dBi5LJ5U83Z/2GIGp7L7nGwzHPNQhS3J7yWaAKe27LkytvA6c/fPn39g4Oqa+fun195VPX3qwLunC2vmH9i/oGZlTdOCgdOm3l0zdZoiv/GASic8yQYLAMhwBiA6Q93NqCLLub9OUmpcstOLaHGCwAsItnQvZqjyadHEUVx6cz+0JMt+sjy645vIQH91edGont0XbPj9msiaPXiIVI2/NHhk35IePbMLh0yeP6V6/ZPPA4KflKlzBqAsnGkVRaCONIPUOstxn/MhJ+nrRKMzxUmcTl2yP92s88eVhKvIfTe2KDHRmKtlyd/2PpPpA3vsPbRzw4w1sz/8snbmA6Or7+w+pUPP8mXDl2wVvqx+wJu//YmVHWb32L5q0oAeXXrkBYa2LZl5056LnkfvwhP6xD0X5YAIN3pyAOvaT85494494cnCD133dnN3O1oEqNZDegiV4IHicLJoMOhs4HS6dC6+LeC2ulLMRKks6LWkMWHX6XqfaELKyMnTOhsGs13PNCxJNkz+Z/0Qg6GhAeewK698pKaNLwyr2caOScrsU1mzMEJygRWCYYcgIoBopDa7TidSq4jaQa/8RJkG7MortqVTEvILI6Z9PL1rzacn//ov0pY1S3t/raYhx5WrKDBA2ED6Yh0dqvitsEECMJuofkCEQsyAJOqq2jzatUOseZR82L1nz+7xMwlZzIVNAOBQIge7xQhgUfrILXa7jtog/71CzQq3qDNoZYbSkOzBpo31obZtOw24a8BDQx4ubWIXRk7UT9S1Kckrtu+bHgSEvqQKP1d3kPleHwFKDSZuX2mGBGlK3sc5EGO7FpnEzw8MXLlQ8pQsvpNv4K4ld9471NP2/hFAoDt1kaPi26q3zgo7lONnEnBvHfMfbr3iP964r4XTTjgzJSYsWHJ0V/3qF3eu3/B8lN07fsKwYRMeGCZM3nHw8LPP7T+w/TH+b/YjjwCBau4hdsY9BF+ZRr1AgMrEoJdu5R/4fBhELEUxdqM72c5aTGef1+IQVnvjPTGxCb3wfhzek01IufGW24c+AOIZzq8gnCYLACAbHrsGKMNHNDV6EPR/osTBA8ziYuCw7Tjs+ThseQz2CwV2Ou3PYeV9xMZBVchkAMkvnuAQM34FFf4CxEZ9KD5qXmxUIBBiM2mNMBxSoY3Sba1zpQWwlbVVwCXk5EIqmmhqKj93lzEgkm2zG3tH7IEWecP9w+9rGZ4ohslCYnXDUm9MGF2J0ihbnJBfkf59Rs7q4vv9Y9X1ozq9+dbRTwPhSMnYbk2zOnXtXqqkXKHH1tZM7NOvw5ip2e0XjzjcWDEhMjB/yIz70jFvcU/eGRvmVKrdoPJ0bltbq9R1v/YaDgTdn4hNzIa84ltA1MLCGETS7SCOQSAGkdoSIv86xGsg3HKMrOsQE6CUQxiaKGmtgtyAkWIwIMNxKIN5QK4xAIk3MIIVnNA/fAdPM+wIOhPaRNEtuvROycm7kHm7iMHM7wabASUqOtByowkglmHm5an5G8bOiYau9y/SAF7vYVQ2zqR5UUeUXdxLDtMT0SMkNXqR9Lhag0cfURpetbZG/AvZr2jRHOZSOkc5ztkqzrMIAf55rM9N5VmbON8PqhxBs8aRmyFqoTwG4b4dxLFrV2MQyS0hsq5DTACHylWC/hhXgUA+gFip9id54Z5wod3t1glmAKcgCUk+rogS11erXC6/JJ+WL8jcIsuyoNfbqiJ6Kri17tNEXW55EDWhHZV7uVhLarxnM5QhVqpNqbM3bcJ9eBf+bn/07S9xNlt4lIyKtaWSunqyntWxHSQcba5nhhhNYrmqS+3jurSmJdWx7jiVLwUx3sKsmLb5bgdRi4YYhP92EMegKQaR3RIiX4PgeGy65RhZ1yEmwMdxnW4b5z7CQrQJJmEDGMEX1st6ino0mXXgy0+0x2rMHLeOu0ewbTh8BHua7RiLw9m2MThS2DCa/3fbaLyfPTsaR+CIsWwrAOXzv877434CJ6RAQFkZnnRvmsAPExtcAA6rqFMCF0+a32f2945YHTpRoDazQHnjnES1lrm3+Fq4+YgL/ygm0lglwc7fxSoM1BZEj3qKzovZ1zsLv1479tEH9ykddGe2jnx04rGmh6Mjpu/9zy/NwbFk68SdWpPhmOUDNr2FDyl9dMMXV699l61D26bmvgOVZjp2ZRN9qTc7xVdOrI9LlUxpXLoVMfk7Nb7fDFELp2MQKbeDOAZzYhAZLSGyrkNMgA3xlRNMtEfCbHWUTvF5CmKjOFSQeO/frHjvH9+pMOtFUbKDBB6vWeALiC8fs96sl2LdkZoVarkRrHVH8v9lCDcaJGexM+zzQ42NZ9GHnuYrO3mL5LvvUdvFy4zXWq/B6ei/V+5Y9yQAqv0oW6R0aK94ppxcMTUAXpMJUu25YkGhw5Hbrl12RaQd5LrV3S5tj+vm0xpaZCBL2vZIQjWCo6Q2/2lnOTKUqE/1UYJv5ZAOKb36Lxv32p+OTCrfUnn27ofnjujZq094yVz2TcPf/v7+58IPi6dX3OnPyC0L3b917LZdPTcF8w/0mVQxcHZN+cTisqHF1YMuXO0r7Nv3562c52pXkOTnPL8TACXovgLUVWlXOH6L57V56vN2t3t+7FP1eajFc/Gz689fe+UW3xc/vP58whegruiOKsCNGRZehzj+cwyiTQwCqAIhKbtXOVDENWdkOJQLre3tedlIaF+WlJTe3ghi5y4pbYNtKyK+AqGgV6RD66BdECyZQU+xzqKriLgsNtBaO9R97viBxZsNL1corarUot3Jy/+qHSkOv7bLFExMz5TiAMaaVIb/wg7NmPnUc0VVb4+a/3xO8a6Hj/0reqcOO967tWbwurHswpy73lz03Mt7Jg1ZtfPpwzvoK7OWGon8BOY/+yddrEUqp/ie+4eMYP/9+yRWGwjyVpav5k5sXH9/5MVNo2XdQ6Sw4ektO5V1zXc4lW4kzreeMU+JFaqnVDtxVIn1ikl8vyqRVppEbn5e21993vp2z4/9rD7PafGcS1R7PsEQk1d7TaLX/gqAo9URXolZHHYXKGOgqI3xIgApTICovZYRgzDHIa79iUMMSoA4xl6IQTg0iG84RDrHQ4OYwA4CqBbHZ9d89VRlx1zyq6euqsJ5fsnUqhXwYN5jsTttkj7YRp9eETFSj91nsfLIR0+9LqSttY3QmLJw6/3b430QyITiIlAqxdlBMcj/lHpUk+6gRVqnV4kwil39+e/sK5T/9sUYXdkp9n3vr4YN77ll3OW+pzc8v7NpC3vppe0vPUtC7Ev2FzR/cQmlWcInr25+cGHXgtrefZ6cNHMlm8b+taaRbXjh4Aku21jXgbraqmOrzaLyJC1RNqNUrt0Vk/1HquySb/e8drD6PPN2z4+p45Ngi+d8fu35a9/f4vtcJtrzCSkx3Wh3fS2Ph2YhR9gJVO1CD4WTPAaDTSACKjsZTifKZjMqJ/QQ8tX1yhOfG8nPjUN6iccXE96Pp8ejezqVFHXsFCrqot3J8iefZP/q3KW8Y1m4nPwYfwOUY3tEGCUsjvv7PvxEa3orl8vQ6iZn76u47uxt1M+b2Kjnf3P2ZWVxBdGcfXw7QXSpTl4Si1SnX6L2X2yaUjNt+Dw0Xd40o6Z25NzmV4rxTJ9pvAljfYjl95r63Iuxboyetf0XbEBQGjL6zuy7cMOvu8aRRcWffLRjTHRO6DzXjNjutSq5e2KSf0PVDI8mmZuf107VNOfWz4851OeBFs+5ZLXnE/yxtZarrfrYDqw6wr2xGWIjpKsAWu+I2t+VyXex0jOkFJfNZpfsrQMOsKeYPHqqT+NdjB7q5euvRZPnb3oYUWsXUUomXo/W9JUVbx7J4HugOKR748Sz333/yd8fMwk63mSElTs38OYRzF9LmyID2Efsvwpjn83sV86KdcDaFQ1NOXQi58u3ce/ZMxo1nF6Nmgn7Y/TmxejV+puEyuv9TaJArLfsb+Iw6gkU6UvxFLggHe4Ot0uSrE5nKpjtqZKY4bc6eDxpBaOR51hGGj+Vwg8UUAc4b5zk4det2ia1fWVJO2TlvZF9aafq7NnSl1EYN4y9zJ7BYRgeN5RaonxdR8+Rfs09fmXXEH+ecs89LqzDiTgeF3ljSZmwlZ1m55QTGn6hNi32qy1yujAU0iAXCmBQuG26zkI8nqx8t7tVlk4oDOW1Mbbh0RHvSCKixdiunWg32pIyxcyKCIieFj7YoVjVRAeseV9R9a0q5rdyvYktTFkxnyvWs/Nzup6pu8B+ROnrBae6djz2+InL0aAOq4Y/e8+QDVf9G154buPm5xvWCb3mrjKRjN+7vp4xEwtQh3q8Y+a0KbPYz19MYDO5tw1mkLIPz3985rOPP/10x9NP7wBEE68Q7pH8YFF6wGWwWXmN0KJs3CSfKkwsE/Igzx1QzhIE0DR3nLfB89CcmUMWLuFF2u+WPJGTu3C+t3TBoiIAgpP5iG2lhdp+kEMyxSpMejflw753u9KSrHUfcfpp29njxj46a8zY3z3YPRTq3rmsqJu4b9TM2lGjps8c3qFLlw78AkQdn+k78TN1N5wPn+Szg2gC/nKrZc73En4mKLYb3o4vKU6BwvQ0olRTQpJEXXkDB/TOLAxZRpmn39tucP/KjIL21tHmqcL5rLZZnbvMquO3Tl1n1aldEci5Ff/FEyCCePMvngykw+K/eMIh5f8VUtYgffQ49lB7+R0HUNTpQenhP6WBBkscHEs5y+QZ1WF29yx63DMUTVyicNM3RdTpRZly061Rq55Od5RisXIk/bGKDPGARzmLjqmfcouq/e4LkcAKAEQZizSpY1khOWwS0KwXbHbQUZP2M1+x3pUgbyrhA/vjeGG9tcNjs9M6maNnb2B4FnXTeR1Tw7TF6DZldL0ZRcHuMIs2WRn9LW10DWe/ei9JQJ4ELUkjOsxJ7m6+QYbnXvbTY2Ow6D6FHh/7lTTBZZSVLOtqB8g4iCCHzeZK+dC1Y38ymWJ3vb5SBnteXszG7cAfyXB6EYzgPBD/URrIP3Wr6u+OqQ9OmDF94qRp5JtZj/9u9sx5C/icym8TiHvgB8gGOwAEwU4c/M4nELJA1RaoJelK5ZPTbBAIlYikk0WuCInpvPM3e2CJ+16ASv2UpGqjUBAIkMRRWhRNSeqtK6QAyGYBkJXxUyYgEkE7ZYLxAQJIVjbPWkkXx4+ZIJRzr1gnnuT0TQ2Xp3rTPZ5kI5Hl5NZ2wZDslYJtjN4kb/+ILklMTUvtHyFp1rT0tPw0qqdJaUlpzsxM6BvJlJ0W3iDhg5ZN3bwwdMsfKruRW2ZQbuRlt9evdcorVpPyolGwuJT/dUDsCHUKOz4AWfRHQvA065Z1snHLxtW7/oddaNewgZANO4LY+n9OPN+rQSxmD80rC7ed1/Rm9/puaEacl3tH9TwUsfXIpYPVzprl6o4iBXdYT0AUtDAtYc3y+EuJtrjkUwGEVlI650ylKvE+5ABA/HNTwuf9lc+BgItUcf0/AgZwQedwuks0ypTyaYjSqY+iqLe60l3E5aIWOZ1mxPuV70toergeGwR4g0v8V2eKi0otVJZJ05xV7GHcsHQO+0ESk9LSjDup6913x/KzVKdeX9THFGzb1v5TDDfpQ45bECoJ9+43cBcf0nCXXr/F8/43notvxJ6rVEnqc1TWG05X9cp+AAQRKWiHl2Knck80KgqljCAC4Aq1QvJpPHP6XaxCImp1FiUv6pwAUXstt2Ud9NrbHGJCAsQx9ufEKktsFtJBzroOMYF9EK/V+GK1mv8PflNJUQAAAAABAAAAARmahXJJOF8PPPUACQgAAAAAAMk1MYsAAAAAyehMTPua/dUJoghiAAAACQACAAAAAAAAeAFjYGRg4Oj9u4KBgXPN71n/qjkXAUVQwU0Ap6sHhAB4AW2SA6wYQRRF786+2d3atm3b9ldQ27atsG6D2mFt2zaC2ra2d/YbSU7u6C3OG7mIowAgGQFlKIBldiXM1CVQQRZiurMEffRtDLVOYqbqhBBSS/ohgnt9rG+ooxYiTOXDMvUBGbnWixwgPUgnUoLMJCOj5n1IP3Oe1ImajzZpD0YOtxzG6rSALoOzOiUm6ps4K8NJPs6vc/4cZ1UBv4u85FoRnHWr4azjkRqYKFej8hP3eqCfDER61uyT44DbBzlkBTwZD8h8/sMabOD3ZmFWkAiUs5f4f2SFNZfv6iTPscW+jOHynEzEcLULuaQbivCdW5SDNcrx50uFYLzFHYotZl1umvNM1tgNWX+V/3gdebi3ThTgVEMWKYci4kHZhxBie3TYx3rHbGr+Pdo7x4dIHTKe5DFn+O/j+W2VnE3ooW6isf0LIUENvZs1gf/LHojJwdpplCP5gn/5gi26FoYa19ZVFOJ6Sxuoz/q2Ti20IKVJdnqvYJwnhfPH/2f6YHoQF30aZaK9J8T026RxH5fA/WPW/8IW4zkpnIfoFLifGB86v0ffm5nbyRs5iaHR3hNBD0HSfTzoPugRM+hdN0x052KoHLBS0tdgpidAiEesDsgWYO73RWQz2LWIwjqnMe/uYISQtlbyf2NlT9Q9PoBcBnrO6I5ELoMeyHkNnIXGdv809H/DXNOTeAEc0jWMJFcQxvFnto/5LjEvHrdbmh2Kji9aPL4839TcKPNAa6mlZUyOmZk6lzbPJ3bo56//Cz+Vaqqrat5rY8x7xnzxl3nvo+27jFnz8c/mI9Nmh2XBdMsilrBitsnD9rI8aiN5DI/jSftC9mIf9pMfIB4kHiI+hWfQY5aPAYYYYYwpcyfpMMX0aZzBWZzDeVygchGXcBlX8ApexWt4HW/gLbzNbnfwLt7DJ/p0TX4+Uucji1hCnY/U+cijVB7D46jzkb3Yh/3kB4gHiYeIT+EZ9JjlY4AhRhhjytxJOkwxfRpncBbncB4XqFzEJVzGFbyCV/EaXscbeAtvs9sdvIv3cjmftWavuWs2mg6byt3ooIsFOyx77Kos2kiWsIK/UVPDOjawiQmO4CgdxnAcJzClz2PVbNKsy2ZzvoncjQ66qE2kNpHaRJawgr9RU8M6NrCJCY6gNpFjOI4TmNIn36TNfGSH5RrssKtyN+59b410iF0sUFO0l2UJtY/8jU9rWMcGNjHBEUypf0z8mm7vZLvZaC/LzdhmV2XBvpBF25IlLJOvEFfRI+NjgCFGGGNK5Rs6Z7Ij/45yNzro4m9Ywzo2sIkJjuBj2ZnvLDdjGxntLLWzLGGZfIW4ih4ZHwMMMcIYUyq1s8xkl97bH0y3JkZyM36j/+58rvTQxwBDjDDGNzyVyX35Ccjd6KCLv2EN69jAJiY4go/lfr05F+Ua7CCzGx10sYA9tiWLxCWs2BfyN+Ia1rGBTUxwBEfpMIbjOIEpfdjHvGaTd9LJb0duRp2S1O1I3Y4sYZl8hbiKHhkfAwwxwhhTKt/QOZPfmY3//Ss3Y5tNpTpL9ZQeGR8DDDHCGN/wbCbdfHO5GbW51OZSm8sSlslXiKvokfExwBAjjDGlUpvLTBY0K5KbiDcT672SbXZY6k7lbnTQxQI1h+1FeZTKY3gcT2KvTWUf9pMZIB4kHiI+xcQzxGfpfA7P4wW8yG4eT/kYYIgRxvgb9TWsYwObmOAITlI/xf7TOIOzOIfzuEDlIi7hMq7gFbyK1/A63sBbeJtvdwfv4j28zyaP8QmVL/imL/ENJ5PJHt3RqtyMbbYlPfQxwBAjjPEN9ZksqkMqN6PuV7bZy7LDtuRudNDFwzx1FI/hcTzJp73Yh/3kB4gHiYeIT+EZ9JjlY4AhRhjjb1TWsI4NbGKCIzjJlCmcxhmcxTmcxwVcxCVcxhW8glfxGl7HG3gLbzPxDt7Fe/gY/+egvq0YCAEoCNa1n+KVyTUl3Q0uIhoe+3DnRfV7nXGOc5zjHOc4xznOcY5znOMc5zjHOc5xjnOc4xznOMc5znGOc5zjHOc4xznOcY5znOMc5zjHOc5xjnOc4xznOMc5znGOc5zjHOc4xznOcY5znOM8XZouTZemS1OAKcAUYAowBZgCTAHm3x31O7p3vNf5c1iXeBkEAQDFcbsJX0IqFBwK7tyEgkPC3R0K7hrXzsIhePPK/7c77jPM1yxSPua0WmuDzNcuNmuLtmq7sbyfsUu7De/xu9fvvvDNfN3ioN9j5pq0ximd1hmd1TmlX7iky7qiq7qmG3pgXYd6pMd6oqd6pud6oZd6pdd6p/f6oI/6pC/KSxvf9F0/1LFl1naRcwwzrAu7AHNarbW6oEu6rCu6qmu6ob9Y7xu+kbfHH1ZopCk25RVrhXKn4LCO6KiOGfvpd+R3is15xXmVWKGRptgaysQKpUwc1hEdVcpEysTI7xTbKHMcKzTSFDtCmVihkab4z0FdI0QQBAEUbRz6XLh3Lc7VcI/WN54IuxXFS97oH58+MBoclE1usbHHW77wlW985wcHHHLEMSecsUuPXMNRqfzib3pcllj5xd+0lSVW5nNIL3nF6389h+Y5NG3Thja0oQ1taEMb2tCGNrQn+QwjrcwxM93gJre4Y89mvsdb3vGeD3zkE5/5wle+8Z0fHHDIEceccMaOX67wNz3747gObCQAQhCKdjlRzBVD5be7rwAmfOMQsUvPLj279OzSYBks49Ibl97In/HCuNDGO+NOW6qlWqqlWqqlWqqlWqqYUkwpphTzifnEfII92IM92IM92IM92IM92IM92I/D4/A4PA6Pw+PwODwOj8M/f7kaaDXQyt7K3mqglcCVwNVAq4FWA60GWglZCVkJWQlZCVkJWQlZDbQyqhpoNdAPh3NAwCAAwwDM+7b2sg8kCjIO4zAO4zAO4zAO4zAO4zAO4zAO4zAO4zAO4zAO4zAO47AO67AO67AO67AO67AO67AO67AO67AO67AO67AO67AO63AO53AO53AO53AO53AO53AO53AO53AO53AO53AO53AO5xCHOMQhDnGIQxziEIc4xCEOcYhDHOIQhzjEIQ5xiEMd6lCHOtShDnWoQx3qUIc61KEOdahDHepQhzrUoQ6/h+P6RpIjiKEoyOPvCARUoK9LctP5ZqXTop7q/6H/0H+4P9yfPz82bdm2Y9ee/T355bS3/divDW9reFtDb4beDL0ZejP0ZujN0JuhN0Nvht4MvRl6M/Rm6M3w1of3PVnJSlaykpWsZCUrWclKVrKSlaxkJStZySpWsYpVrGIVq1jFKlaxilWsYhWrWMUqVrGa1axmNatZzWpWs5rVrGY1q1nNalazmtWsYQ1rWMMa1rCGNaxhDWtYwxrWsIY1rGENa1nLWtaylrWsZS1rWcta1rKWtaxlLWtZyzrWsY51rGMd61jHOtaxjnWsYx3rWMc61rEeTf1o6kdTP/84rpMqCKAYhmH8Cfy2JjuLCPiYPDH1Y+rH1I+pH1M/pn5M/Zh6FEZhFEZhFEZhFEZhFEZhFFZhFVZhFVZhFVZhFVZhFVbhFE7hFE7hFE7hFE7hFE7hFCKgCChPHQFlc7I52ZxsTgQUAUVAEVAEFAFFQBFQBBQBRUARUAQUAUVAEVAEFAFFQBFQti5bl63L1mXrsnXZuggoAoqAIqAIKAKKgCKgCCgCioAioAgoAoqAIqAIKAKKgCKgCCgCyt5GQBFQBPTlwD7OEIaBKAxSOrmJVZa2TsJcwJ6r0/+9sBOGnTDshOF+DndyXG7k7vfh9+n35fft978Thp2wKuqqqKtarmq58cYbb7zzzjvvfPDBBx988sknn3zxxRdfPHnyVPip8FPhp8JPhZ8KP78czLdxBDAMAMFc/bdAk4AERoMS5CpQOW82uWyPHexkJzvZyU52spOd7GQnu9jFLnaxi13sYhe72MVudrOb3exmN7vZzW52s8EGG2ywwQYbbLDBBnvZy172spe97GUve9nLJptssskmm2yyySabbLHFFltsscUWW2yxxX6+7P+rH/qtf6+2Z3u2Z3u2Z3u2Z3u2Z3s+O66jKoYBGASA/iUFeLO2tqfgvhIgVkOshvj/8f/jF8VqiL8dqyG+d4klllhiiSWWWGKJJY444ogjjjjiiCOO+Pua0gPv7paRAHgBLcEDFOsGAADAurFtJw/bt23btm3btm3btm3btq27UCik/1sq1CH0I9wl/DTSONInsjxyKcpGc0VrRNtGx0dXRF/FpFiV2KbYl3j++Jz4vkTaxKjEgcSXpJzMm6yb3ALkAnoCV0ARLAcOBjdCAJQJqgWNhJZDT2EbbgTPhz8h+ZFJyDbkFSqgVdGh6Br0BhbFFCwHVhNrj43DXuH58V74WcIkahHvyDRkLXIGeY18SxWl+lMHaIVuSc+h3zHpmNbMJOYuy7DF2E7sFvYMJ3Clf+3DHecNvjm/m38g1BYmioxYS5wqbhZ3S0Wl2tJkab50U04pl5CHy9vlmwqlZFJaK4uVnco55YlaUK2kNla7qEPV6epi9aMW01jN0zJohbRZ2mptj3ZWu6e91wE9vT5LX63v0c/q9/UPRiZjprHS2GmcNG4ar8yIOcycZC4yN5mHzMvmE/OrhVq6NcCaYC2wNlgHrAvWQ/t/e6w9115r77XP2fecrE4xp65zwM3lNnZnuBfdZ17E071sXj6vrTfP2+Hd8F74lJ/eL+Hv86/6D/23Qfogf1A+qB10CAYGk4LFwdaf2C+JfQAAAAABAAAA3QCKABYAVgAFAAIAEAAvAFwAAAEOAPgAAwABeAFljgNuBEAUhr/ajBr3AHVY27btds0L7MH3Wysz897PZIAO7mihqbWLJoahiJvpl+Wxc4HRIm6tyrQxwkMRtzNIooj7uSDDMRE+Cdk859Ud50z+TZKAPMaqyjsm+HDGzI37GlqiNTu/tj7E00x5rrBBXDWMWdUJdMrtUveHhCfCHJOeNB4m9CK+d91PWZgY37oBfov/iTvjKgfsss4mR5w7x5kxPZUFNtEoQ3gBbMEDjJYBAADQ9/3nu2zbtm3b5p9t17JdQ7Zt21zmvGXXvJrZe0LA37Cw/3lDEBISIVKUaDFixYmXIJHEkkgqmeRSSCmV1NJIK530Msgok8yyyCqb7HLIKZfc8sgrn/wKKKiwIooqprgSSiqltDLKKqe8CiqqpLIqqqqmuhpqqqW2Ouqqp74GGmqksSaaaqa5FlpqpbU22mqnvQ466qSzLrrqprs9NpthprNWeWeWReZba6ctQYR5QaTplvvhp4VWm+Oyt75bZ5fffvljk71uum6fHnpaopfbervhlvfCHnngof36+Gappx57oq+PPpurv34GGGSgwTYYYpihhhthlJFGG+ODscYbZ4JJJjphoykmm2qaT7445ZkDDnrujRcOOeyY46444qirZtvtnPPOBFG+BtFBTBAbxAXxQYJC7rvjrnv/xpJXmpPDXpqXaWDg6MKZX5ZaVJycX5TK4lpalA8SdnMyMITSRjxp+aVFxaUFqUWZ+UVQQWMobcKUlgYAHQ14sAAAeAFNSzVaxFAQfhP9tprgntWkeR2PGvd1GRwqaiyhxd1bTpGXbm/BPdAbrFaMzy+T75H4YoxiYFN0UaWoDWhP2IGtZtNuNJMW0fS8E3XHLHJEiga66lFTq0cNtR5dXhLRpSbXJTpJB5U00XSrgOqEGqjqwvxA9GsekiJBw2KIekUPdQCSJZAQ86hE8QMVxDoqhgKMQDDaZ6csYH9Msxic9YIOVXgLK2XO01WzXkrLSGFTwp10yq05WdyQxp1ktLG5FgK8rF8/P7PpkbQcLa/J2Mh6Wu42D2sk7GXT657H+Y7nH/NW+Nzz+f9ov/07DXE7QQYAAA==) format("woff")}@font-face{font-family:"Open Sans";font-style:normal;font-weight:700;src:local("Open Sans Bold"),local("OpenSans-Bold"),url(data:application/font-woff;base64,d09GRgABAAAAAFIkABIAAAAAjFQAAQABAAAAAAAAAAAAAAAAAAAAAAAAAABHREVGAAABlAAAABYAAAAWABAA3UdQT1MAAAGsAAAADAAAAAwAFQAKR1NVQgAAAbgAAABZAAAAdN3O3ptPUy8yAAACFAAAAGAAAABgonWhGGNtYXAAAAJ0AAAAmAAAAMyvDbOdY3Z0IAAAAwwAAABdAAAAqhMtGpRmcGdtAAADbAAABKQAAAfgu3OkdWdhc3AAAAgQAAAADAAAAAwACAAbZ2x5ZgAACBwAADiOAABYHAyUF61oZWFkAABArAAAADYAAAA29+HHDmhoZWEAAEDkAAAAHwAAACQOKQeIaG10eAAAQQQAAAICAAADbOuUTaVrZXJuAABDCAAAChcAAB6Qo+uk42xvY2EAAE0gAAABugAAAbyyH8b/bWF4cAAATtwAAAAgAAAAIAJoAh9uYW1lAABO/AAAALcAAAFcGJAzWHBvc3QAAE+0AAABhgAAAiiYDmoRcHJlcAAAUTwAAADnAAAA+MgJ/GsAAQAAAAwAAAAAAAAAAgABAAAA3AABAAAAAQAAAAoACgAKAAB4AR3HNcJBAQDA8d+rLzDatEXOrqDd4S2ayUX1beTyDwEyyrqCbXrY+xPD8ylAsF0tUn/4nlj89Z9A7+tETl5RXdNNZGDm+vXYXWjgLDRzEhoLBAYv0/0NHAAAAAADBQ8CvAAFAAgFmgUzAAABHwWaBTMAAAPRAGYB/AgCAgsIBgMFBAICBOAAAu9AACBbAAAAKAAAAAAxQVNDACAAIP/9Bh/+FACECI0CWCAAAZ8AAAAABF4FtgAAACAAA3gBY2BgYGRgBmIGBh4GFoYDQFqHQYGBBcjzYPBkqGM4zXCe4T+jIWMw0zGmW0x3FEQUpBTkFJQU1BSsFFwUShTWKAn9/w/UpQBU7cWwgOEMwwWg6iCoamEFCQUZsGpLhOr/jxn6/z/6f5CB9//e/z3/c/7++vv877MHGx6sfbDmwcoHyx5MedD9IOGByr39QHeRAABARzfieAFjE2EQZ2Bg3QYkS1m3sZ5lQAEscUDxagaG/29APAT5TwRIgnSJ/pny//W//v8P/u0Bigj9C2MgC3BAqKcM3xgZGLUZLjNsYmQCsoGY4S3DfYZNDAyMIQAKyCHTAAAAeAGNVEd320YQ3oUaqwO66gUpi6wpN9K9V4QEYCquKnxvoTRA7VE5+ZLemEvKyvkvA+tC+eRj6m9Iv0VH5+rMLEiml1XhzPdNn3n0rj6/EKn2/NzszO1bN29cv/bcdOtqGPjNxrPelcuXLl44f+7smdOnjh09crhe279vqrpXPuM+PbmzYj+2rVws5HMT42OjIxZnNQE8DmCkKiphIgOZtOo1EUx2/HotkGEMIhGAH6NTstUykExAxAKmEqSGMFl6aLn6J0svs/SGltwWF9lFSiEFfO1L0eMLMwrlT30ZCdgy8g2S0cMoZVRcFz1MVVStCCB8raOD2Md4abHQlM2VQr3G0kIRxSJKsF/eSfn+y9wI1v7gfGqxXBmDUKdBsgy3Z1TgO64b1WvTsE36hmJNExLGmzBhQoo1Kp2ti7T2QN/t2WwxPlRalsvJCwpGEvTVI4HWH0HlEByQPhx468dJ7HwFatIP4BBFvTY7zHPtt5Qcxqq2FPohw3bk1s9/RJI+Ml61HzISwWoCn1UuPSfEWWsdShHqWCe9R91FKWyp01JJ3wlw3Oy2Ao74/XUHwrsR2HGHn4/6rYez12DHzPMKrGooOgki+HtFumcdtzK0uf1PNMOxwDhN2HVpDOs9jy2iAt0ZlemCLTr3mHfkUARWTMyDAbOrTUx3wAzdY+niaOaUhtHq9LIMcOLrCXQXQSSv0GKkDdt+cVypt1fEuSORsRUwgrZrAsamYJy8fu+Ad0Mu2iYFhexjy9FIVLaLcxLDUJxABnH/97XOJAYQOOjWoewQ5hV4Pgpe0t9YkB49gh5JjAtb880y4Yi8AztlY7hdKitYm1PGpe8GO5vA4qW+FxwJfMosAk2X9n9X2cVVfnA36pzHNHJGbbITj75NTwpn4wQ7ySKfAu9u4kVOBVotr8LTsbMMIl4VynHBizBEJNVKBAfMNA9867j0InNX8+ranLw2s6DOmqIHBIbDfQR/CiOVk4XBY4VcNSeU5YxEaGgjIEIUZOMi/oeJag4mEB3PUOweCaG4wwbWWAYcEMGKn9mR/segY3R6zdYg2jipGKfZctzINQ/vxkJa9BOjR44W0OpTKAskcnjLTcKyuU/SVIWSKzKSHQHebYW9mfGYjfSHYfbT3+v877XhsIwGzEUaleEwITyE2u/0q0Yfqq0/0dMDWuicvDanKbjsB2RY+TQwOnfvbMUhiNPFyDCRwhZhdjE69Ty6FjoOoeX0spZz6qKxxu+ed523KNd2do1fm2/Ua6nFGqnkH8+kHv94bkFt2oyJj+fVPYtbzbgRpXuRU5uCMc+gFqEIGkWQQpFmUckZe2fTY6xr2FEDGH2px5nBcgOMs6WelWF2lmiKEiFjITOaMd7AehSxXIZ1DWZeymhkXmHMy3l5r2SVLSflBN1D5D5nLM/ZRomXuZOi16yBe7yb5j0ns+iihRdlFbd/S91eUBslhm7mPyZq0MNzmezgspUUgVimQ3kn6ug48mntu3E1+MuBy8u4JnkZCxkvQUGuNKAoG4RfIfxKho8TPoEnyndzdO/i7m8Dpwt4XrnSBvH45462t2hTEX4Bafun+q8jIzK/AAEAAgAIAAr//wAPeAF8egd8lFXW9zn3PmX6PNMnPZNJMRRDMkzmDYgZMRRDCEmMMUPJIgZEepHlRYyIiNhRUdYuS4ksy9reLDYsdOmLLC/Ly7L2CgKrrCJkLt+9T2YyYPl+D8804J5zT/n/zznPBQKbACSTvAEoqJAdtUhUJpQYjBJVAUrKSkIOJ1ZUOEKOUGkfV8ARiPB7E72m87WJZF58ibzhXPVE6QsAAnMufI4H9XXsUBh1UpOJSJLmQNWqNsasLkKhsrKnA/T1HCF9PQzSAPYtD5V5PW4lmFeIK86EcCRbObLp2lGjGxpH4+f0wLkjjU3NDSNGxYSMxbSdDkzomhE1SypQalCISvniob1lDuTL7injC1O+Mr/xmeJtxeRt/iJviJ8mmrjFOr0BJCZ3QAbkQFu0ypCZ45HcRqNJQkiT/LKsOO02s2Ryudze7CxVUnw+v9+tmKTcgEEymzPRlgN2e5rHaeOXyeeiisnJFagMOSsqSkr45kL8Tr450SfM5/y1V66pGvBwTV1BcYcDEX67QjQkbo8cigTplyVI2OHh/6zdXHO4+iR6SjoxMPzo8O21h2tPx7O2lmylNV/tY5Nwubj3fXUA/8BuFveBr74CoNB84V6pSnFCLhRCL7g7OijfR7Oy3FalR49AcXYRFBnsQUcgkAYO6H15j6wiAGu+I+Ao6pleFDAWKJZMX+aImNunWOpiskIVH796ewAqEzvV9gqX9nQ4Qd8S/1V/ScSM/rmsTP9FfNUNIvzuVlRPMFxY5PB6fY6iwsJw3/JIOOTx+lT+WzaR+xYWecrR7fWFFanqi/33nnn9+v+MvXr7mk933/v5Gy3PrN6yZjg7WFV1D5s2oGoh7nx+k2vvTrkeDT0HKlieXvvakkfecj/5uKnhm6iNHRk27a6bevTL+clH3ulVkX3cBTJUXjip/CDvBiO4wQ95PB6qo/len0+WTRpofo8nLa04mB3UgpeX5PbMLEzzKz4/tapOlXt5a1llpXhN7FF7r8zJ37o/iN15Q2XhvsE8RdajOqwFyrwFGETXr/0F9u9dNnZsWW9869X1azow9qe/kpc7D52mPRf//HcJFrR1npvf9sWX336EO7/9x7lqeUMn6frt8y+//ZD/JjzecOGEAnxvWdzjpTAzWtHbGjRhlhdMXqvLVZSWnl5kpSoChLJVtcwXSPea8vNLSrT0dEnTegyPaZIUqIlJLnSKhAV/pfBuhb9EbE53bYVIM/3S45hfiZ+7th8IFPHN5QuXcscms1vF8kiAZ2qBsEEEFQX7FnJDeNy+8nIF2JLZ7/77DPtk3rJhVV9vefPD+57CzCF98cr82+s631s4/vbxrKPf1XjT0Iqrh/+uafTMxR+9e++mxqZnxzzx5l8embstxo7PeX0Ju3DjoqYJA7C611hyd3hAtH/zpD5jAAVm4DM6Zjj5C5WIAIu9DuxCIB0kuvEBAKGBbSTz+L+3Qm7UZjaZqCSBqtrN+VQgmAMTua3joeaMhBTicTt9wULS8PSj5x58eNk9Z5c9RUrRiPte3MTKzvyHRd5Yh9vFygP4yq3JlfmyfHG+so1LyP/5yqgRNVjuDPclRSGvk7Q+/ejZJY89/OA5sTT7ifVb+zru/OEM7tv0EisFhErSJGUpbrBBOOo3ms0ypVZUVc0umUyqilarYrDxpN1aJrKQuykJwvwz/yPMUOCTXSqlRa6CiEzJy8U4J8DWf/jpM/eeOMZeLMKpxYqbPTyx088Oz8MKtnMuFqefm4gzAKEZPpUqpG1g5qivGRSjkSKAxWo2giJRKOFCysqS4vjNhQXCAa4Bxz1HEI+yNlx0FBextqOk9SjezW49yhaIHbGzuBtOggKe1wgFWVapDCXbdSNt5ghfoNCgMxLA3X1v++dV+eg/vIsdR9MJYWVcS5rISqDg+CuVQQLkSiTc7QoHPANIGq49dw6wi7GwgmvujZoUrrSRNsaMLqjsmfjnkYu4aU6SlJZ28xECNyqt0mMrM2pBricBidueiNS5iDcRA0ir4h+y4yQgGJP/DwLVF05IQ+W9XLoPLou6LYoTFPCnGT0jYkaV2kfEaBok8y+1kkYCeeDQnIEyQI2nUrlDE3kkDT3PzsfZhXMoxZHGw2OmTRl7w+SpLeQoW8gexttwNi7C6ewO9hD7/usTaELr8eOAMA+A1nJtTNAj6jJKAAZEs8WgqihJRgX9wJHOkYoXkf8iwR2RiKKqRRiitWw3lYdnr30cDzNae/8Tw/1L3sS5gFALINXpKDQgmp1pQxW86M3O8aoqMTlNtTGnSjATM2tjXEgCYfS3hKyuCkFHkzBeScI6WKhFVxLuD+EQLt4TkOo6CU5f1drrhvrrVly/dspDayfe+8EtQx7fuJG0HcbZLyyc1r+5qXbojtE1xa0dt4x/5c31r9hA6MYtP5DrVgijoiV5Po6KKs3MBOCVStFlgez8bG57v8/vq4tZ/Gilfr8pX7VqJm1EzJQGeg3j5/xX8ruWMbrG4oduFyXxMEFyQlkpkMeJTvhKbCMY1j/o2ykPlEmSr335KxvYPvbZydev29P65KNrX58+c92zfxv6+Kil76PnU1Sl6fe+l694//zIweMjUO1ZPnH2TU3fxqa09+l/6OHXAQgEAaSZuhddMDiaZ1epkRAzpTKAxyVzrnGh7JLreGi7qF1VqO5WvoGQ0DwF584uo3cpz4sCBzc9T9SAQPKgoqI082X2QfxhshCzXmZ5Jmoo6MvOYAk7gCWH6cudN5+98oSroZZNBoRWbuEw1ygDmqI9OZ36aJrbbTPYqIFmZrldRpdFA27ONADF4/HXxjyKYhkRU9LgYsIJ6e+pgHAkGUjkgUhLSBg2N9w3IMwpylMaKScT/n6efcC+PLN8xActmMGOhu+4bH6EpsV/yAgOoO0n9/+HnR2B5h7hr455LAPJ1+wc+1i1AYGhXOs6eQf4IR+uigYUp8WSlweZTnAWFNpz6mJ2u4d60kbEPGnUwENEvUTbVJbqTCjIAQJlPo8IXEUNdQEJcCAhMvd/gvy8Q3E6TmsbErv++Z2tRuuN/7f1X+zsNyv/vYhoN066sbVlcRuZiq/iWvuP7rEb/7LuhyPfsFPLMffdxfMnz7+1fu5qEc0RPdM6QIHLo14FgCDKRFYNMiWU1MaoAsLfupYpQwobhpDby4OfkoJ4iZQWPyy9jNLm8wLSdEtUyzvBB3lwOVwbLXYqnl6U+o3+Qo/Hnp1ttBtL+ihOZyBQXGwBS0Z9zJIGwfoYXGwTYYlLnVeWdKFwoCSqAj0/LqoW8qk7kShFiku3kK9cfCPVHyDedt/qpeyLL06zk4uXtU1DyfXfE2fPmrng0Ccjbhg+flxtq7zz3ZUzXhrU/O6sjqN73mrbXD2iY/Kzm89vbBp7Y/3VcwaOI3vqq674XdnlYysH1Ym8GajvcgekQQFURnOzZJfFEgyCCwqLtNy6mKZRrzd9RMyrUkMdR+Nfdbfu7DIBzCIaw0J5kS16edcXuNOdBXwbyU1J1ewxtvTOqxtHP/3+JIOl3xOz3v0nmr9Y+f2d8VNjp4xrbbm7jQ5mdazJdtYzasufW2r+83/H0fEE+3DTXbdNum1+Hfd4stOSZuvMURh1OXnyAPjtnsaYXeumMPAnaOwXTOb4NVYT72PqU+xG7xcf6mPNQAQX6/IUcHKmcllV1UUlBRXFZdIaYyZNUjgzJ6Rpm8u6mKrApzM0vUgYbrTrbF2SFHbS18Xa5GhSmF5P7JYqZODSiqKajIK/VYNEqQIEZRigFxShVFwJURhGD6JU0ZlDP443kvW7ccNSPH2abWFfCns140peoYDeNeZHHSqlRgkMcp00ViJSV30QKhkjagSue7JMQH4304/FkrTgKC9Tjh69VLueUScBrhFPNVAUJJTKEur6Ce0u1dCFuorNZH28UayJb2IaDjjNtKWsWmioXPicrpB365FYFc3LTU9PA+B2dlqdhUV2QCMFCAazGmNBl900ImaXkg7mVCR4KJVkyfpRJFR5F86oRckaXOFoe0m/7W6YevPVY5uWvzf1w3P7vm99YGyIHU4139VjH6ob1tLvqqpxR9u2r5m2onVI9RVXsHUX9eMTLkxQdnCc6AuVEIv2VCsq3G5XOGzt77rMZaWBtEDvNOgN0au8hkhEMg3QTPzqkVUq5feAklS7rOucMleiPU7ivc6kQtuiYCqrfNTdlVF8fxLxCKgtj3iUQC44+jrzOa06UfyDSESH3x2j106vnpWmTXnhlT1o+UfT/qt9NdGau79/Zhf73+exCP2T2Pz/ZefZXez6I/gIyv/EkRs7Yf3IFpM1FG27n5x++NQ9Q/otPPTGQSQBH/Pd/9Yf/vjjne1sx152gh0p6f3eKHwYW3/EZZ93sA627uCCpcfMzwj7AIC8WN4IKljh6miAWKkBQZHNZgqip6CSZLOSmpjVSs0yBZocIpTouZRiZWGortKL8gsDiITjI5Uik+LHJ7FXiYTziRJnywoMgWdwNFstbzxXRcbikdvy72CqiPvXAaQznI/t4Idczsm9VLdbktKzzeY83vfZ7QGDlqalDY9ZNLRSTbODPb0mZneCvyYG9BLcSxY9KQVDSTe5ArmSp7voCQYwWfE4HPqnwOu4AyOYNn/C/fPZh2fjx7C84/aZ8xev2nXHraxT3vDKpkVrHaacdQ++/xGdXTuy8Zr4NrZo3PgNgDCXI/UBnh9eKI36VZeLN+NWnxscUBNzSKpskmtiJleyNBOvSfVEKuQRD2+0Iw4l2BUdoTI+ZiikBS+9h9OfOtrxL7aJvdiOkQOHDrc2tEs72U/HmW846xyGi3DSZ3j9azd1FvUDImwoz+E2NIBd1OtGAIdVkjTZUhOTqWTlLbMzaamUcEELnGVzAbVA0BHKleew8ew2Ng534wR8gL3Dxq5ZjO/xGuQP7A55A7ubrcHDnUMBdY8RLs0Mg6L5BgnAqphMiBbFWBOzKNxLAnII3zehaKqJofOXXkp5iCsitPAkbol0bqDV8RN4ijmIm4tl7zK2BLqkUsalGqFvNN1AqVkBQDQJoSl5QlZS0MVSLhaCX7P9dHD8OHKMEwKWxLu8KBdxL6ZDTbQo3e8nNquVEFemy2DIsGlmjQdbOr9BNkt+r+zlsmTu1FB3wd0z5VlnstgW8BBwKLpv9YJL5RlPdMKNOALkU1L14E93sr+yVfg43vTxgZtW/GXnd1vevKGVHafhuOnyAlyMU3AcPjDybB377rOT591Y2mUHeYJu/Ug004jIzW+QJFm2GGhNrMaABoNsUijK3QmbMnfKFN2XPIHtjr/NdmE5uRrDZG78Xj5t2EIGAOCFiawBT+ozgRw+bSAGXiPLwM0MRsr79e4NCw4Rxa5IJL6kRnJurq0bOKEZy79hDV4k7gVL5JHn1l4AdgYS+tfxVS0wMJpjIcRkNiOAzUBl2cq/UrNZoXwP3VtwpgBXF1eWAOXEQAdVfSMRDKBcx1awhYvEZm7FB7CZETKxJf4D39CN6/Hf8XkJ6VIlly6LPUkqBVCQArccJKJUl6GXoPq6r3PD1MsbzldfSPxvRcyR3dAvmukGo9nI1bbxUPHKisdJjEQxq9QGilBcN36X0mUp6hA6Y9DpEYujXuXykscVRBpkK4wudhzbcaSC07GdfUgtRrZEms9Wzok3cw1WSi3nqklH6R3oPr8kYcedOm6WR9NMYETFagVwUFlRVM1MVW5RVLtHv11adI/EnAKwL1KEcM/JO9nv43fpSiwh81U7+qQGdrQtXseFv4FZvycdQPQ8+VKfDHgE0jgAfBZF8RpdNTGjRO01Mer6daQROSBexQQy16Hxpkj+kj3BXubXE3gz1vNr/PlDb76Bs9nSNzaSY+xxdivejVP5tZCj0mP/OYvf4smfoAvtpHU62rkEFkhGowdsNrvdbQXBV3ZNM9TENGr/TSzoRn/ZLXHoEyAo4ckJSx+au+BBspEdYacX8yA6iCb0UGXmlKkTd504Fz8rb/gchAXYat0CdkjjEZynUFmSCDVIJg9AhmYypVOVEwBXRFK5UWSV22N7Ev4uHU92T9OQe+LX7PPaKziWzWZnfL9pJMZW1bO5OPS3LSUP1S3lg9poocvnk0ySppm8njQw8cTzu4wWMA6PAZgtFm40C/WaRcikzJbSWfPzuXKqQ0sxKLdfgl3BF0A82brsgaXLW7gB12EPzH7oTqxuZWvZKtp73M0Tm+Pz4vvlDUeOLdxZwVwPk1KRVS2cQX0ce4s4n+RlpKcHICC7LeCGy4rdAbAELNlGX3ZNzCdRYyq+uhvwVHHWrRpn+IvGGoVFl/MhDadWMcJP9LZen9cr+din7JuOx/ZeN2FqnzFL7767DtWvZu2f2TrnyermlsJrn977BC7f/lkz5g4srx3e8+orqypveeqmzf8qL/13n8KGgcUDKqrHbRP6FwNIYiqrimdLCgBFNBhVKlHOuxSdv3y2lARgcoLtYrOlOn53IGEMEF7k+dXC13JCQdThQHSbDQaX08hRhsdSYuuXVBAOtyLx4BHI6+6CYLnlEXbyLfYFex/D9zz7BAf0ztqVZ+7EwHn6YufCPz33/DraBqjXfyHBI2K+RonRKAOiVZYkC3BDJ+q9VNpUJOaj+sXtVx6h57CC2dmLTMMKdPlKFXO0a4DY+dTwvZeN/qJLhrqRy8gSsx+T0e52yQh+v2ynlszMrKwci9mcnemSzdRvt6NJiOSi+EtCbgo1UyM3WkiKOMKJUtMlGvCIi78nPihD2fPbzWFJ6WPdxqngfix9q9Sr9HQdwoJDth5mUy/nm1hKoRixV/mpUJxwVT85trLi1EAa6twb+aS+9uuhNBsStmnSbVMVzTXLnPpUo6oYTYpJ0C2VLGYDkWXJqFCUkhDL9evG+ooUZ3VpjZj8Izex59h6fnXg56wfNmF/DGMtC5Pi+GHyHdka/47Y4j27dJCYyF2B7wZVlZEQEERvNFFF4QqiSgVDdslOjEH5Z65AarLLowIDZAGWchEZbA/LwDo6mozsXBTfQUqoXleVJiZ0RugfzTJISFUVEExmlYuSRP1I0IAGUcZdOgxNpl1qFqqPbALSzPPvkbfjTVJ6vIrs30m/RXi/0ykkLWUbyWw9T7KjVgXRIIFRJlTBfN2EuvH0BNZX4iUpmc0y8bOPPmIblXMHz60Xa1gA6MDkVFt/ZIKYnGpfnBa6sUmAHY9/mJhqI4S4fJ+QL55xoKIY+VYNoOZTiaaCvQtCfCFHMMy1CH34IX7GMmfKjQd/UoR8AzFIA+R3QIHeUTdBWVYkSTznFd6SVJko0DW+xLKLeyTRZYcwiGjADQ/jqVO8uP6KGOiGzmqyKN4maq1OtpHWXhja9SRIRonoRhEaJZ5K0NrOFyl//vMAAGKNdIQ+qATAwK1gBjVKRVTIdwCUpB/rioP0XWLww7EvHPD6PGRL5ZkqbKpcLx3ptW2gZ/z7GYIdmjju9pfm6E8Zq6OFTovBQvLy/P78LIMhaEkbFrNYZLfbPjjm5jWdnDM4JnvBk0Az/y+ZVYSeXlcUJWdMvMcN9+1u8h0omny9N6YT+huGr1r0xzd+Or/5xbv/On7T8Y9PswO/X3znY5MWPHHDsNfXvfono1K6rn7f+K3vx32E27h55MJbxwOBFVznDsUNTsjh7BvIojRg1Mw2n89szrWA2WPUFFDSh8QUL7iGxEC7mCz83SHi7H5mUeZ0aISzRVANCgTlw1AfH9d2D8WobftHX+7YNsMT+hpLLZbJM2ZOJJNvaZk+Q5rNdrPv2XH2t6XzFTdbPuiJ9jP3rwh0PPOXNWvWAMLoCyfoMWk2eDi6esRYymclxCubh8RkDexcM++lZZJuOTk32SdwmnJoYkjgUBQyIf4DZqJx81Mjh9525cmTzcuHVf/BTQZgFvauOZFVwBH49ZIydr4kH4iQK81M2CcaDRi9Gi+obTZhqFy7xwIOIyi6fTTdPt5ft4+oT4Q+ecShOXlPGioU/BLkji3iOnVPiAnZ9vHnOw9ON/mw7Jv+1omT5kyVp7dNmDnLjWVoRx7zq9vG4YSfTjyy5vt7ViWNk9BynD61y+DMEKROSUpzOLKcJlOm3+OkzuoYFVUUVMesmuoZHFNTel5aloiry3bI3RbgrbNeR4XKwOMJ6AVAxMMtOP2GaQZcT2aVs+/Y3zDt7LdoiJfID985vmNc3Qb61PyZM+d3NmAPdGAahth3Jx+789Eel5+4rCjB7nSOkgMeuCKa7SZElSn1+qwAPhndyHVz283akJgZqJ4bgp8v7QVDiRwWFgxH9KfOeieocBWpiZ1l+9eu3bj/ufm1o2uv6ocGOq9zCZ23rKHh3ZdLPsoafsVgoKAwtzSV26sYyiEKd0SrzFlZAwZIfRwOUqzmSkGUpIHpPXr4fJFg8Kp0K1jRqlj7qv2GxYy5Eke5wr7FpDpWXFxYWDksVqi5e1fH3BkXz+n4pxIOWz79gRHv0LneqJs2FQ76ewKfPao+pSsqEvmsj+ykQFfCF6ZeRcGFyUQK8v26El/4WGzqS33OfxjpXbL2ndc3sTfYvm9+vP3WksHVg5tvOnmsZKGTFc2buvrNabOfa5w5/drrmura10otT/ceNqZjJ5Xzew187smt/1i1bPw9We5Roeh1xYVrZ732vkM6L1UOHVlb2WcEHT5q0qRRuwBhBYC0lmeDB8LRdATw2Y0Wg8Fo9Nolp1MaEnNqJkCjR6D/JfU5336yUOPaKqJJEuCQeFQirWX7O+6YxfZjqapqE/61bQ958LsXt8S/40CwpeDekav/vh0ILAPAD7lsA1jEZFcyGsFksprtJg9Rr4kR6DJ/ZWoO7uobKtNnnyJUlrW3X3ttO14phMgLHn98yIjzPqkFgFxoY259XSt4oSTqd/L0JgaDT/NcE9PAaBctOk/sjOTEKYEwCRGJxwB6tajQpMDBcxoHXzN8CJbum6GLZe60066mRmnd+eJXN6mThXRIWPMH/Un+NdGgxLmTUKrIsmYzWa0Gg8lkN4P41WCzUcXkofbu2oTf3cjSZdpuokXRuGOyi1dx22KswGZWhYd5AffOIrF9jYxdh40sI74Et93MVivueDXr0gYPcG0ouF4DRIkAevQioLvExgPivyvuhO7qQJ5BQRgeLXS7XPrsKDMzI6PAajSaTPkuq9WRKzu46XwOzWzPRJNH7+G7krl7+OC8ePqbjJDCRIiEfKFykdziVfBd8q+ke9n++uvnTGL7vy529F437Xwso/dL097ZwvbVXz9jOnlw3rz12+LfSS1Lh1+/urZpy+F4kfhtxYuQjGCut1tMFxHAq6vrscoOoatQFU0Xx29SyV/XLRG8TS0ierkyof+ZtWWXEPbn7boC9dce3JHE5yf0pzhpostXLJYMcLnSvcYhMa9mp0Nidu8vu/xUrvPeVQMOCCQs6MzrxGVT5986ecr8W6dQmX3ELvzxh7swGyl/I6Xt6/70Qnv7mhfYKbbnQTS8jE7s8wA7B4LrOep1cC1ckMMn1Hl+RVFNlKpZmqrlcuQEq9U9hBOEwa5mQEaKzBKmSBWoSQVlTvPepDFCnPndRKFJtuemosq2GZrG9p/taZv8wfaPbt58TGf7vePdSx/wsv5K9SPtbB87/T/s7H10mU722JDgM67pTN1euaIq8dIsyh+TpOUZ+fg6PcNnz/ZanE5V4I0FhsQsv8m6iSfIBUmS5S2dL8HBXl8ook+LIkFBaLdMkafPPzxZ2v7R5zsmPXeFIQMJ22e1lq48uri9oOMZ9uLa9lNYiho3Z9+6xqU/bcBDAybXN3ZFFJ3LddVEh0mcejw5BCxZZVnUS7wGFxqlMrTMRy+JIqpdWewrCD+6iu3/sre97yvSbCP7xLR8SXyH1LKxZTYkqp/1XIZ4dpmjpLktAEU5bnchWNw5lhxTli9rcMynUdPgGPX+vJ2/2BgiqPTHK2HB5clePsGgXCkPt082oetPnbx1/bDrDtW395oycuG8yJd/3/Xu6MZHa5Zcv2zRrf2wZn1HILfzsvKx+b0rCstHz73+8VXN/8y//JriK/qHR/+30LeE6xuRa8AjToRYDHa7y2UyEIfB4fWZnHbn4JjVYrfL3HVyQt3QpktOVnRhgnBcxKOXvoLpIyFPwCO6cjK3bsas9tdeeHRt8xasYDuu+TD4aeiNN0jGwgknTn4e//yqK4UOT/Gc4zM+cENZ1E8cDrfby3t/j9NoJ7JNtumyPcmJ1sVDgItr7tQYgH+grxdrpR2zt72PpSLjsXRp7XUHt5Mj8dki4Ynt/EpI9JkPcrlm6BV1m0GWiYgIK0G0GNEuC5llKWndDU1X/x0SbTfiOtaElf/INyryZYexkjVJLfFF86aMXUzaumS4AZRtXEaWOMsoSyaOIVng81ETVTMyMjNzVEXJ9plMVLbbMxQ7yDqidR3RdPz2LIDSIO1WQ8wBsin/pGskRZpuUfew19lm7LMwJ1eRcrT7sG6R5NCsqBgvN92NPdk7uARPdt4vtTDH4m9q1lxH/PGvvE03jMkcer4XnuKKI5gApOW6bWqi+YoMaKSUSAQlGWWzQVWtfIZmMSoUAA1mj4T2S2cBqaROkYZeq3KlhdkClOu/mD2BI48cxZHsMWxja46fYO2kPwmyZ7A1fiy+DRewhcJLzK17ycs1KTC73ZrXK0koahm/Jgob/pNT8no0p9XJMTHDAFyVskQJkKKvhBlTUzxHyokifvTqgNsSaw9mmBRz7n4cwoqu+vcfR9RErqqfl+fkfr2/YcZNo8ic866XXnR8Z72xNZI450HXce2MIn+oKqkIYDYgmvQhAm8c7YR/MwyOoefSIULSSMJGySlCWEwR6LrOB4nC0uhAZiCmDrLp6+3xekDI4T38Id7D54ipCHUbcnIcfn+uNTMzIFGXy8qjKd9qSbTzYosp2hbbF7bnuBrm+REWRw08Coc18VTQ4xFQ6+EJhDmL2m6/c/OZG4cpn31T3XpmM9quH32qucGAVz7Z9jEdXMUObcyzBF8xskNVg+knbU8BIO5gJWSlYgMK7tcIpZJMAaCyhONDYlbqCOKOo0cV29lA1ylOauB7yBN7yOHlOmgGQ75bkoI52TabW3Z7qCzl/3/2IIuHzuFynuSi2BZnlftyiBSnzxyCyzwcrImh4e0Xbhz2+9mfKtWtL7xTP39x26LeM2aFPyFVQ7CnuWmyw5K3EXsOrqIfh2dPY5tNjY2nGm7QTxGQIqmCtoEHIlG/Ag4zmKnd7qNeu82mSJSaHQ5QoCRU1lYi9ElBdqqp5pwa1sv/RAMmELwQB0baym968pqFwxaOC99ePv7pgf89chFZcXX5l1NzcyPRii+nphf8lzhBwpbiQanl0rP6Dg26zurbad4v56mukCugE0Wi7Vh7JsTasSV5lIO0dJbKBcljHAhLOdJqfN6cwad7QYchPV3OyCA+n4mYMrPSXCNiBtuIGMiGNH4pGWmKygXqpwH4S8+ePzvOII575nOCTh4R15lS69q26gmSEBt94OCr7YtF6z7vlm8b7mpdcN+rL/fHcyhjZk77c8arjmflv/Bn9kZObzbAuFFEB4A0ST+d2BztZXeaidFqTfd6iV/zO51ado7Fn+avjxnT0sDFqcleG3P6QR7xs+NNXUfUIJTSVqjbjT+pBpRfbpXXFSKawsFwiBuQbNyyZcyzs2sbcS679w9k3/mvbhr+6qufy7sbvojGrt10dOm6WtZ5ttes1keObtl5BAjMBCYFpHXcnkW8R87TLC6j7EsnBrDZ8jIhM/OyYp9LSycWo2xQPZ4ctYBHz/YyHc11H2qb9S+iA4oURXyC3SM+0WGqPrVIoJJaFCmMXFRdbixfuGzBqEk3j1qwfGE43Pbogt+Nn93Y9siC8v1T6+qnzxxRO50cnPC7BcsWhCMLly6MTZs8uu2RtlBo/iNtYyYOnz6ttm7aDBHpCoDEp+PghZnR/7I53U6Plce2UaYyMYkJqxeRED/HBp/idDkbYkCRuuwmm93WEFPtdgt6FMsl5xX9mtiW3kNfypcpEhAfkgPKkCfoEXdAGF7cGCBD0YAVbOGWH374gX38448/vsOW4BViZBv3vHrfq8eO8RdyHMhFiKNCMGoniiKGmUaJSlTVsUcEbCpFdAhyJGBIAFHnAbag8wAAgUm89lnw/0o5D7g2jvTvPzOzu9KCJNSFaAKEBMYHAokSuQpiY04OODjYsWxCcjbkNaluuPdyiXuaS0jHpPfeE0N68fVO/ObSe+8uy39mVlqEzr76oeyi+bG7U3bK83yfkUZBGZwCMyKlaRaXRRTLC6E4JyfkAld4DKmpsbkrK0ttpSafxzc15nHqTVNjepQycUvmivi5NiuyMYtA0qyNo3NOVr9OFfZJmt75WUW7VMhOWtE4fsubj9zRP33SzuaW6LxFB3rWTJj4xSuvXdHyYsOAb/bpj257c+OS5s4tvmrim7appHXPputbn8kPlVdURssit194/xklXdGr7p3261Hh7uKKUGH0uu2nzi8Pxya1V5qmAUYu4UfygiRwVi0/YrQaWIvIdGcQ4pBB7dzU9snCdpLZJF/SOXJNjdRPPa0uMhVd2TKurqk5Mq5FXFPXEB0/7ucNExvqGieOb6wDIIw7lSbR99oBPqhmvm9ikm0mm7/c7yzPc+bV1IrpYEmnX1mlhbZglpActKMVbEo36zBrHWyifBGnSASrw44ZvIhr6bwgFCxiuH4R45HIul+c91p4c3j55tf/fvilPddGFx5b8zJqf5X9DCi9v/m10vvcrj6U09uHsg/0Ke/29invHSBfX7VJ+TAv99nwkcNvfNd82xjlI/4/Su+rLyi3/ObXaPaLTJb0b6xlBfCX+DHKMLqgAOoieZk65HLlmXXU56PLK/RmGI2e9HQbys4GEGweShSEA0F1mAtak3BQbR1SPGxVVo3K6irbp3YM1ToJV3pGr452r7n58XnrWi6tr79h3tY9yqTy/KbYvMvxsYvGRLrPu/BCWegef0l+cNcmpeGP/qIz6oqkNPas06Fd6BEEkMAIbZHRaUaDTKd2RMKCgERqGDdkGNkrBpBGCE4XBIMoIpOMsR4lWko4kLBqJI+K5j8Faab66Q897w8yR4ALIR3yqYfpaPGg8hFyDSo70RG06A12/oayC49HL1E/s9K3DL2QNXzKGb8fhTCZCCJkRZgzSkcQkogAAdYJoQTf6LXQWZQQHjx2hLz1I7pgEIaGErEHWAIzAAhaezTEW+S5kUqBYFHUgcViJEbamxB9uT/ROLFE8QLBIegdsp5+naSN8spKbara53ErgY4FlFnoIwadmhP5X7VaYcvuz5QHAu8h/cO3K+s89eFTJuceP+dft9utd0xUFqDpyj3kqh3K1+H6uhrlzX/ZctHQEckuSNLhJG8MjPTGCNLRbwWDZH+Fr/6Jm7D5hAmyIDMiQ0ZGTrbVkMkqRQ3FUq17vL06HSowmDyctbXd2N5201ln3XjW5a88G6uvnz2nLjJHWMg+7W0766bZL10emd02YWJ7G+NFAYSwiCGdcx+ZGTqdRB35BoSomd9sMRrSZYQkAYOKeoYC8S5MM5WnxriwyfZwnAs9I2/h3kG0RVlFY12UNylYiiCAo/gZTriVRKwOA5LAgiyuTNnkwQ4Hyucer4lJXb96j39EPHUF+JnjK/5+briipGXeqiuf3np9+4YudA6O3jbYEQv6S2bt37Cle8be7rMBwVgcxo+Ir4APJkRy7enY7QbIl/LTzVK65C8mdrvDIed4PSa5IIE5pbQ8dlABTRX6S6xu1DgHrezj3QjuuaN9/n1P7N541ards5oXtJ3REgwFWsOdE/b9v3W9wlu7a432i6at2N7wzOzzq6tvrAr76ePuDExYn+qLI0JEDyCnCdwXdyjui3uFjR/VNMjMIUk6ao6YiGZWHZ0i/DX75U5H1aEgAOK2LmrkhkxmMUmXJFnOsjrBQR/drXNlOGl7yiCq4Y2Z+zTTkbYwT8qwtv73xo0CxS6XhZtDZ7WvpVaAD0ZnlC6fNWF+vigy+yj67YoVdz/PrAF7Z8wo/9mM65SDUhQQLFSOCbslO2RAIOJINwsiAoTMFr0emUykKWYSWc8XiHtk4gMlbe5qgAb7UsMIa0IFwu6bbumd0PqX1/72IW5Tjkmn/3QfCVmPHEWCwiKd8Cj0e7KGEUURmUU6Ebk1RiCQCHSypSLhfEr/+2Eqe2hQsaNeALBCVcRlNjI7Fh1Y7Gaz0W60ySYW9pXNXt9QQI0EXB1/3PjAIiZPQYprQ3RWgnr3Xd88KXuOu/GW5v7s6Kwj6xc5btOZJpzh7hmf2cktXDiKGxPRSYI8MjopD+WfMDoJeePRSb4QbvyciNkVzReismdxFD2z4Oyi0vHr6MwOwnTUfEt8ic9KPBFjIvYqgzhkDw/xTGK3kxc9YlKPgt969IarH3/wwP4nFG9dY+PEiY2NdULbnf0v3Hr7wAu3dHR2dnTMm5cy6s2OlKZTy49OL2AW1Ib01FNiGh70BD7YIdHEB79/Oej1B9UBL+6NL0aoFonqQehRdg4ip/LxIFqsSMPn2KuMXYbaUNsyJZw1fMrGrnIA6Qpa2n5Y+TuAYvg1fgUA6eAP5Nrjj4L8IMFW+uJUVye0D51Au5h8T7W6B7CZSZlyNlXeJ75ClUs8XEnM8as+Eb9qmXpVwDBeWUH+LLTzNU5DpKiQug4YJk0jh0pMoyDbnI1lQp0JPk9rzJdhoRy8xZvKwaN4g9Cm5HHsnddbrUub3bCVWHLF4ldiF1wYPjM27aFzzp37w3lvHP3F7rOrUcnw6jY6d1dT86yJ4eiY0sOnTO6//YLru+j0cyyamXhHhoZU2lu3GPuhiOexHiQ0HfQPYqfoh9HVJ1B0w2//heIgzFQV2SMV52iKgYTCOlIxU1N0cUXaQwR7uWRYkxbXSNDfPYvXhpfEa4MpdD7OPtrg4sg4yUbMNmIRLCjNZEJsvgbgEETRbiYUvqb4syENGQkj/JFkkzkxTAQrMmlscsKiQLvUAAeUNb8G7yQ062PCs0QKkEYsI9rR6nzH9imOvcoLeLew9/ghbKIUT+hoLlq5jiPvcYqZDnXNrC6WKXZGjNP8+VlGYAXOBfY556p5+ZaodTT0KC89ZE+UXqqiG9pSFPdShT1JcXDoO1XhHnmNmZqia+gnXgMYFag1wGbucZ7cAJnQGCmivUCW3ep0GlBamtthAIqVWwGovcRJi9eKLYy8TgmP0+BgddahWmkscQqUlpiPo4MhBwPPA1tV5FzFz7cKwm9+d+CzzzahATIdd1Du/G5GoOPWnR9+ofQoyl1qHsRXeDuriLez36eUA+dUeTlUxtt7N1fgvJMpulHDv1AchOdUhXek4hxNMZBQZI1UzNQUXVzB2vvoeGkj2IAMglnogXTIjaRLBGTZYORGZXcgqMUn8260FqnLBlSM7lL+uB+Vocqr6Rhetkf5tfL7vfj3qKxH+SMavZf++VuaSiUAhD7DLeIHkgA2yIZCCEdyXJ4cuz0tB9LAW+TMK3Ab3QxXJQWpdOWImbyK8arGGFaJqpEG2V2IO/yqihEFV1Wm94Xts3tnv8iA1RevaL1x1sDRP56CjrR2UWL1/ZBiOG0+WqzyvXWXXHDpANrEwNWGNfM3DSi/fHYJ/rbsp+8e6j5uKR4aUmlIXgO18Vocrdaz1uOkKrqR6V8oDkKPqsgfqZipKbq4gr0RJcl9kqDwq4yNv3kb1KtYuCSJSmbrqZpIDiOjjbIoSpJTMDbFZEdTTJAFWdIRyZowKGrdjOZBjePIDroW0tZGwh2UUz1yNcPaH1CQ4fikjst3rbt0NcHv/agMUij5c2Vc18rz5/NZJM3JfMkD1dAaGU3tegXFxQDlWSZTbXkgUGPKKtBBcbEui2SWhkqnxEIQcFgyozFLwnGq7ZUx0g03TH/aTYLqcnOkuuX8iaFL8zhXsVAn4a3SSDRSWl1/RVfoo3fmXTau+ubIbfnTo2vnNjQ0TVjXsWQjbb4+hL9FfuGvkV+cNqai1JldVTJn7srmu+7JLfy6KLhqVGhcaeOylsh5lbWnl49r6TrnKPVMv/LO/azH5ASbVEBr5VQ+UtQfAPb2jbbEazY1vfvCE6Xna+kHfxhi6RUj001a+kAasPTikemClt4lAX+3T+GCYcUDmqJ/lKrwqwogTCEpQjeUQBBOgS2RydU1JDM/P2g3GoNBuabG7/GMKZPlsC/fW50fjVVXsyDp7OxQNJZtNo6aSoF3p+S0NFDHPHgbYiBJgQZGv/ERLZmZ0t5q6wkJKnqMhzBz8MufZG0ZXsZRzHYYrWJk1TDShwoZfiVWbn2rce4L19/03NdfPRtr2nHzvKc/emdx/d3LDyM4XkaJq+cfm/bY8bqFq1fv6FyOvX+1oHvwefbOru7Y0zcz5q91cn3Tq52bInXKZx9RCGvWp8UlOEsQzpxD6T/05acLVrNap952xtZhP0xWx0+0iY+fnCrjtT1FbQ2389oqStRWanr34n+eflDP00eNTBe09C6rWpeVidoeugYAvcGv8LTaXynTgF0DGRLXuBwA/y5J0T00eaRi6JdU8UmS4qDyuqqwJBTvUMXlkqApuriC9Vdu9UkSBIfk5fPVpZGx4MYuV46oJ+kEY0tOTnr6qEKLpcQNmZh+SJ2ImdjppB56CnnSKS02+RpiJifBU2MEnYC8izsQ2clwI9I+1YYLf3Gtkw8SVgdtm4XAwyNdtX46hDAvXCL2GCmnN3ZetuitjjuuvUr5/0PfKX9DwuFDDfpT17zfga0rz19x8fIFq84TXdXF99Wdtr1n/m5lz4fKh8pLyPrJR8gyV+hdtuva4/Mv2Lj1ih27+lg74MwMf2tPV9/aEPAZUHI97ucl3KK2k5t4PReeOJ319ZfAyRW8pRiS+gUt3aSlD6jpeSPTBS29y6C2pIDWK8yCw0JYeIl7wbKhNGJ1pqWZBQEIyYUcNwVKAXHz0vPBYdBQiw8WTxJRTWOGj2+K1tf/PFpXNzVaf2ojO+KOwcEvTpva/POG6c1EmNrUMqWhpRkIfcaHKAN0OZ81eEfOGnzxWQOjb0jBFAZx/C+zhmCNsJ9hQWsvOLVn0n5GBm1eUrt/zK5jR21o/OiJKy9AhwzKa/6alefjSoYJlXV2dVyL7IwUqpp+Qes1ytH2RjTouvnWlnFKMOP2oSGVpeD1c2ZST4ByefGmpvMavgVOruA1XMnTC0emC1p6V0B9A0u1np977PkV5qi9zXh+BQ8XJOgmziYWsLhqD+1vHQZzli2Dxi8VWsCcbXDIRM6dEpOdxEnL+CQocxLLTDtnDWdWTT4Wyh0nAU7ot8Herhf//uZLf5xv0ulUfvGjOONEDrXMYEgzK+CtE9qVsXpQVixvbB7mnLQ8CVqeut5Qc/0zNdcJKk9oH6byMk5M5VGJGk2mO108BE7wQmekxuJwGFF+vs6WAeDL0umKLHa6drMgI7HQX0YznaWSNBddcwhCLotpRQ5tBcd+ThplmiAy+BMMx2M6XcOLuERnVGvx+3WnH9vn31Wm9Cv3oTPQhPGbvaRDW9Q9dstdd/XVrfR7t8jpaBvqQuejTSZZXeCR145+8+1PDivZbnPyN+hT3SphMXhgNARhQWRMoMKEHQ6/X19RkWu3V+Xr9aEchzvgiMYCATCbfxaNmc3YJNDOmfLEZnDT4VwQvFNiQupwHj45Cp00iOdT56kG4bniI7dDo6KTeT2fSk+Ltyhf7dl5pPfHLSgb4QUvT7nsi2+R+bhTt2fL+U90tDx99FwN5Pu4fbWMBnC3/ZprdiD9/ciByqY1XcvYaf26naXlbOCeHGf7BhavuJhFHD0h/FXwSAVgZP0Zi5ozAMh6jE0ZWF4vsh39sg5pyx2NKqQzEZ2XGU+dFNAgrdc1Ne977elTUafn6kbhr2ed0XJ29tMLqh5sYBENqFX4M4lKD8Q9ehmS1eqmkUWyR8ay7CDxvRTYHVKNZ7qk8YhEdy1YcOklCy+67Pqa0tKaiorSGvGlCzavv+iCDZu7ykKhsrKqKkDwa+HPgkEygQuqIm4KNEUEQjLdBhvobPTrYvM6MzavFyCQ9fpZmoNENQebXw6qkISXvbF5mNVHiE23yjF6xRM27knfvXTUtKZoET+/fAk7F+uray7vKyjOr+KHAr4bGHqI3IN7+G5S+AS7SU0nbeih999Xlbp/qtQllG7Sj/p4jIw7kiaIOqTTySBou5KZB5gLq7jGWhvCumKTs7N6sN5L+p1zkG2h8t3HkHQFCVwRmQhIknSCRC8wvD8WUrffQHtNwbWDkz3iI84XlPdRySFI3luLeVIwEfnuWhIEtNuffHstwOzeZBl/+gzwRczUIGsiggSSZNFlkHRtI0Z+oT8E+bOoWSnwxY/oUzVPdILhSZyRP8ezp2Vz+E4SGJn/ndpNDXwrMFMaMYjsRi+qN9Luoz60qB5QH885cqO31JNM8Ua1DBJFgVlJkOt5SRihMGIaeQcIpN7Ap91gROGgt0eWkkvbi2wunXrfKIyCdLA9wszuRplAgHssUq3uc6/avnXvvku37cGf9hzou3r/LbcAELbTizQXhfm75mXsYF6m6kEvys4gbKuXAofMQuS5LUhtbJnmP9AJy8gdX3yp56m7v+Aps89kZzPacGPqPmctKUf+VkA7vpHbtCsijrgDV9RLQAg9pa0JI9VZmsxW0W/VN5vqlE12xKZeO24nRzp2bfoHPRPEf7z2SBs4vvHEBm8ApCxj83oe25YVSSeAEcaCFtqW8B8j5EX48mN//IKMjge2AeK7BW0S+6EYdkQaJaL3+XI8RW5ntmywWIrSafaLika5cnP12dklBpdLzpRy83Knx0heRt66PJxOMvMy82yFPiiEabFCndlkMzXHbNp2YiNNoxZenyxzKUghO/CtQOhvro/H5DgKdA420DrVfS4oWELdb/7qWvq7BuL7XXhXXu9CVyrtGKN5yj0hZNq9ecn93ynPj9q6VMBLtvjQpG+e6ps7ebnwys5f3ucNFDzwTXgIxqK0Tx5wFVff9zVyT//Q4+XsWgfzjp+0n6MTYDbdHRriMbs/Sh7wQyNfQ04lboD45x8nfd7MPgcMBhzF34tPQRpYGbthFXUmWnBEBixim90k62TJikTRaiW6PJLPDTwBLSYu4RpNwn+8DhpfWI1CfA+zWrZnHP5+zefKBrTh0zXKHkmuzliH39q3rwfXHT/UN3Nu1gWuZ9Wn05u0pyuGRuJWn14KAMTT4QTpzcPp0q6k3PF0dS8BvtMDAcsjIIiIQGKXQLYPAt8FgTU2uvZ8EQDruB3sL/EV7krVDmZIWNNupYoPkxTdQ3NGKoYYgS4mKQ4q76sKS0JxHADfqZupKbq4gq9wuaT6/wCVeR0IAAAAAQAAAAEZmiehT9dfDzz1AAkIAAAAAADJQhegAAAAAMnoSqH7DP2oCo0IjQABAAkAAgAAAAAAAHgBY2BkYODo/buCgYGr9zfPv0quXqAIKrgJAJZXBsIAeAFtkQOsGEEQhv/bnd272rZtG0Ft27ZtW1G9dYMiamrbZlgrqN17M89K8uVfTna/oRs4AwCUGVBCU0zQl7DAlEIZWoPOfhXUs0BbVQAL1CG0ZepQd9STPdUW9dQ61FGN+U5LpOW1pswUpmU0hZj+TGOmWnQ2lPNyV2rEoO/A+mUw0CwATG8cNjkwyXzEYZrG9Of5NUyy+XBY7Q4Hm9a8tgCH/WU4bOcwPfmsjc7GvDcYPWk7StjU2G8qAf5xwHQE6D+zHRXUbqzi96bmrEQNEeim4V965jWnB+ho0sNRHnTn7E5H0V3nQAlaAGsawqkxWKfGhDPoO2Ts/Gdwsk5fIecd011vh9O/OaegHO9toBWAfYLM5JBSxvoNquliyEeDvUucbeXvMd55vIqRtTGMJTnzAkP5bdnsXvTX6VGOPkbfYe+yRgh/6xHoLms6QDmmlvyFPThTB2PEtbczfMbr3XUu1JD7fmqUjaYre68jzpPD3wJIH6QH0RyQ5L6Ui/GeGFqDOZLiPj7iXnpkDsKJ5+TwO3LmEe8JYecb2fcazoXMC/Ed4z0J7EFS3MdH3EuPJJX07gom+ff4/DMcpS1ee85bBLQNGO84cgiqPerpVcghUBEeK/S1jzBBfUZbwUv5X/7bkOlslqCEwJ5TBw4lBFsBJdRuHA4vYk/own8RLYvLrQAAeAEc0jWMJFcQxvFnto/5LjEvHrdbmh2Kji9aPL4839TcKPNAa6mlZUyOmZk6lzbPJ3bo56//Cz+Vaqqrat5rY8x7xnzxl3nvo+27jFnz8c/mI9Nmh2XBdMsilrBitsnD9rI8aiN5DI/jSftC9mIf9pMfIB4kHiI+hWfQY5aPAYYYYYwpcyfpMMX0aZzBWZzDeVygchGXcBlX8ApexWt4HW/gLbzNbnfwLt7DJ/p0TX4+Uucji1hCnY/U+cijVB7D46jzkb3Yh/3kB4gHiYeIT+EZ9JjlY4AhRhhjytxJOkwxfRpncBbncB4XqFzEJVzGFbyCV/EaXscbeAtvs9sdvIv3cjmftWavuWs2mg6byt3ooIsFOyx77Kos2kiWsIK/UVPDOjawiQmO4CgdxnAcJzClz2PVbNKsy2ZzvoncjQ66qE2kNpHaRJawgr9RU8M6NrCJCY6gNpFjOI4TmNIn36TNfGSH5RrssKtyN+59b410iF0sUFO0l2UJtY/8jU9rWMcGNjHBEUypf0z8mm7vZLvZaC/LzdhmV2XBvpBF25IlLJOvEFfRI+NjgCFGGGNK5Rs6Z7Ij/45yNzro4m9Ywzo2sIkJjuBj2ZnvLDdjGxntLLWzLGGZfIW4ih4ZHwMMMcIYUyq1s8xkl97bH0y3JkZyM36j/+58rvTQxwBDjDDGNzyVyX35Ccjd6KCLv2EN69jAJiY4go/lfr05F+Ua7CCzGx10sYA9tiWLxCWs2BfyN+Ia1rGBTUxwBEfpMIbjOIEpfdjHvGaTd9LJb0duRp2S1O1I3Y4sYZl8hbiKHhkfAwwxwhhTKt/QOZPfmY3//Ss3Y5tNpTpL9ZQeGR8DDDHCGN/wbCbdfHO5GbW51OZSm8sSlslXiKvokfExwBAjjDGlUpvLTBY0K5KbiDcT672SbXZY6k7lbnTQxQI1h+1FeZTKY3gcT2KvTWUf9pMZIB4kHiI+xcQzxGfpfA7P4wW8yG4eT/kYYIgRxvgb9TWsYwObmOAITlI/xf7TOIOzOIfzuEDlIi7hMq7gFbyK1/A63sBbeJtvdwfv4j28zyaP8QmVL/imL/ENJ5PJHt3RqtyMbbYlPfQxwBAjjPEN9ZksqkMqN6PuV7bZy7LDtuRudNDFwzx1FI/hcTzJp73Yh/3kB4gHiYeIT+EZ9JjlY4AhRhjjb1TWsI4NbGKCIzjJlCmcxhmcxTmcxwVcxCVcxhW8glfxGl7HG3gLbzPxDt7Fe/gY/+egvq0YCAEoCNa1n+KVyTUl3Q0uIhoe+3DnRfV7nXGOc5zjHOc4xznOcY5znOMc5zjHOc5xjnOc4xznOMc5znGOc5zjHOc4xznOcY5znOMc5zjHOc5xjnOc4xznOMc5znGOc5zjHOc4xznOcY5znOM8XZouTZemS1OAKcAUYAowBZgCTAHm3x31O7p3vNf5c1iXeBkEAQDFcbsJX0IqFBwK7tyEgkPC3R0K7hrXzsIhePPK/7c77jPM1yxSPua0WmuDzNcuNmuLtmq7sbyfsUu7De/xu9fvvvDNfN3ioN9j5pq0ximd1hmd1TmlX7iky7qiq7qmG3pgXYd6pMd6oqd6pud6oZd6pdd6p/f6oI/6pC/KSxvf9F0/1LFl1naRcwwzrAu7AHNarbW6oEu6rCu6qmu6ob9Y7xu+kbfHH1ZopCk25RVrhXKn4LCO6KiOGfvpd+R3is15xXmVWKGRptgaysQKpUwc1hEdVcpEysTI7xTbKHMcKzTSFDtCmVihkab4z0FdI0QQBAEUbRz6XLh3Lc7VcI/WN54IuxXFS97oH58+MBoclE1usbHHW77wlW985wcHHHLEMSecsUuPXMNRqfzib3pcllj5xd+0lSVW5nNIL3nF6389h+Y5NG3Thja0oQ1taEMb2tCGNrQn+QwjrcwxM93gJre4Y89mvsdb3vGeD3zkE5/5wle+8Z0fHHDIEceccMaOX67wNz3747gObCQAQhCKdjlRzBVD5be7rwAmfOMQsUvPLj279OzSYBks49Ibl97In/HCuNDGO+NOW6qlWqqlWqqlWqqlWqqYUkwpphTzifnEfII92IM92IM92IM92IM92IM92I/D4/A4PA6Pw+PwODwOj8M/f7kaaDXQyt7K3mqglcCVwNVAq4FWA60GWglZCVkJWQlZCVkJWQlZDbQyqhpoNdAPh3NAwCAAwwDM+7b2sg8kCjIO4zAO4zAO4zAO4zAO4zAO4zAO4zAO4zAO4zAO4zAO47AO67AO67AO67AO67AO67AO67AO67AO67AO67AO67AO63AO53AO53AO53AO53AO53AO53AO53AO53AO53AO53AO5xCHOMQhDnGIQxziEIc4xCEOcYhDHOIQhzjEIQ5xiEMd6lCHOtShDnWoQx3qUIc61KEOdahDHepQhzrUoQ6/h+P6RpIjiKEoyOPvCARUoK9LctP5ZqXTop7q/6H/0H+4P9yfPz82bdm2Y9ee/T355bS3/divDW9reFtDb4beDL0ZejP0ZujN0JuhN0Nvht4MvRl6M/Rm6M3w1of3PVnJSlaykpWsZCUrWclKVrKSlaxkJStZySpWsYpVrGIVq1jFKlaxilWsYhWrWMUqVrGa1axmNatZzWpWs5rVrGY1q1nNalazmtWsYQ1rWMMa1rCGNaxhDWtYwxrWsIY1rGENa1nLWtaylrWsZS1rWcta1rKWtaxlLWtZyzrWsY51rGMd61jHOtaxjnWsYx3rWMc61rEeTf1o6kdTP/84rpMqCKAYhmH8Cfy2JjuLCPiYPDH1Y+rH1I+pH1M/pn5M/Zh6FEZhFEZhFEZhFEZhFEZhFFZhFVZhFVZhFVZhFVZhFVbhFE7hFE7hFE7hFE7hFE7hFCKgCChPHQFlc7I52ZxsTgQUAUVAEVAEFAFFQBFQBBQBRUARUAQUAUVAEVAEFAFFQBFQti5bl63L1mXrsnXZuggoAoqAIqAIKAKKgCKgCCgCioAioAgoAoqAIqAIKAKKgCKgCCgCyt5GQBFQBPTlwD7OEIaBKAxSOrmJVZa2TsJcwJ6r0/+9sBOGnTDshOF+DndyXG7k7vfh9+n35fft978Thp2wKuqqqKtarmq58cYbb7zzzjvvfPDBBx988sknn3zxxRdfPHnyVPip8FPhp8JPhZ8KP78czLdxBDAMAMFc/bdAk4AERoMS5CpQOW82uWyPHexkJzvZyU52spOd7GQnu9jFLnaxi13sYhe72MVudrOb3exmN7vZzW52s8EGG2ywwQYbbLDBBnvZy172spe97GUve9nLJptssskmm2yyySabbLHFFltsscUWW2yxxX6+7P+rH/qtf6+2Z3u2Z3u2Z3u2Z3u2Z3s+O66jKoYBGASA/iUFeLO2tqfgvhIgVkOshvj/8f/jF8VqiL8dqyG+d4klllhiiSWWWGKJJY444ogjjjjiiCOO+Pua0gPv7paRAHgBLcEDlNxQAADArI3Ydv7Vtm3btm3btm3btm3bD7VvBoIgLXVVqCf0ztXT9dzd3j3cvcX90CN5Snmae/p45np2e356gbeH94HP8Q3x3feH/X38NwJwoHigQ2Ba4GBQCK4NfgxVDE0OnQr7w1nCI8P7wi8jdqR4ZGzkRDQSLRmdH/0UqxTrEVsbux/PHe8b3xh/lgglzESJRJfE6MS6ZChZJzkj+RouCA9GJKQuMhI5hsZRHR2A7kZ/YZWxldhtPDPeFd+IPybyE0OIy2SIrEy2IneSX8mvFKB6UpfodPQYeiOTjmnK3GOzsCPYpexaLjdXiRvBHeJ+8BX5Lvxe/qOACmWEnsJ60SsyYjqxiLhE3CoeE6+LL8RvUlRqJXWThkszpJXSbjkq83JaOZ9cXm4gd5IXKZACK4qSSSmiVFWmq0lVUtOr+dXyagO1oxbRSM3UsmnFtOpaC62nNkqbo7M60HPppfXaemu9j77X4IwUI49RxqhrtDWOGzeM92Y985lFWWWtcdZia4d10/piU3YZu6+91j7rME5xp5szGVAgDcgBioDhYDpYDjaDE+AmeAW+p8R/A5ajfCcAAAABAAAA3QCKABYAWAAFAAIAEAAvAFwAAAEAAQsAAwABeAF9jgNuRAEYhL/aDGoc4DluVNtug5pr8xh7jj3jTpK18pszwBDP9NHTP0IPs1DOexlmtpz3sc9iOe9nmddyPsA8+XI+qI1COZ/kliIXhPkiyDo3vCnG2CaEn0+2lH+gmfIvotowZa3769ULZST4K+cujqTb/j36S4w/QmgDF0tWvalemNWLX+KSMBvYkhQSLG2FZR+afmERIsqPpn7+yvxjfMlsTjlihz3OuZE38bTtlAAa/TAFAHgBbMEDjJYBAADQ9/3nu2zbtm3b5p9t17JdQ7Zt21zmvGXXvJrZe0LA37Cw/3lDEBISIVKUaDFixYmXIJHEkkgqmeRSSCmV1NJIK530Msgok8yyyCqb7HLIKZfc8sgrn/wKKKiwIooqprgSSiqltDLKKqe8CiqqpLIqqqqmuhpqqqW2Ouqqp74GGmqksSaaaqa5FlpqpbU22mqnvQ466qSzLrrqprs9NpthprNWeWeWReZba6ctQYR5QaTplvvhp4VWm+Oyt75bZ5fffvljk71uum6fHnpaopfbervhlvfCHnngof36+Gappx57oq+PPpurv34GGGSgwTYYYpihhhthlJFGG+ODscYbZ4JJJjphoykmm2qaT7445ZkDDnrujRcOOeyY46444qirZtvtnPPOBFG+BtFBTBAbxAXxQYJC7rvjrnv/xpJXmpPDXpqXaWDg6MKZX5ZaVJycX5TK4lpalA8SdnMyMITSRjxp+aVFxaUFqUWZ+UVQQWMobcKUlgYAHQ14sAAAeAFFSzVCLEEQ7fpjH113V1ybGPd1KRyiibEhxt1vsj3ZngE9AIfgBmMR5fVk8qElsRjHOHAYW+Qwyumxct4bKxXkWDEvx7JjdszQNAZcekzi9Zho8oV8NCbnIT/fEXNRJwqmlaemnQMbN8E1OE7Mzb/P/8xzKZrEMA2hl3rQATa0Uxs2bN+2f8M2AEpwj5yQBvklvJ3AqRcEaMKrWq/19eWakl7NsZbyJoNblqlZc7KywcRbRnBjc00FeF6/enoi05EcG62tsXhkPcdk87BHVC+ZXleUPrOsUHaUI2tb4y/8OwbsTEAJAA==) format("woff")}*{box-sizing:border-box}body{padding:0;margin:0;font-family:"Open Sans","Helvetica Neue",Helvetica,Arial,sans-serif;font-size:16px;line-height:1.5;color:#606c71}a{color:#1e6bb8;text-decoration:none}a:hover{text-decoration:underline}.page-header{color:#fff;text-align:center;background-color:#159957;background-image:linear-gradient(120deg,#155799,#159957);padding:1.5rem 2rem}.page-header :last-child{margin-bottom:.5rem}@media screen and (max-width:42em){.page-header{padding:1rem 1rem}}.project-name{margin-top:0;margin-bottom:.1rem;font-size:2rem}@media screen and (max-width:42em){.project-name{font-size:1.75rem}}.project-tagline{margin-bottom:2rem;font-weight:400;opacity:.7;font-size:1.5rem}@media screen and (max-width:42em){.project-tagline{font-size:1.2rem}}.project-author,.project-date{font-weight:400;opacity:.7;font-size:1.2rem}@media screen and (max-width:42em){.project-author,.project-date{font-size:1rem}}.main-content,.toc{max-width:64rem;padding:2rem 4rem;margin:0 auto;font-size:1.1rem}.toc{padding-bottom:0}.toc .toc-box{padding:1.5rem;background-color:#f3f6fa;border:solid 1px #dce6f0;border-radius:.3rem}.toc .toc-box .toc-title{margin:0 0 .5rem;text-align:center}.toc .toc-box>ul{margin:0;padding-left:1.5rem}@media screen and (min-width:42em) and (max-width:64em){.toc{padding:2rem 2rem 0}}@media screen and (max-width:42em){.toc{padding:2rem 1rem 0;font-size:1rem}}.main-content :first-child{margin-top:0}@media screen and (min-width:42em) and (max-width:64em){.main-content{padding:2rem}}@media screen and (max-width:42em){.main-content{padding:2rem 1rem;font-size:1rem}}.main-content img{max-width:100%}.main-content h1,.main-content h2,.main-content h3,.main-content h4,.main-content h5,.main-content h6{margin-top:2rem;margin-bottom:1rem;font-weight:400;color:#159957}.main-content p{margin-bottom:1em}.main-content code{padding:2px 4px;font-family:Consolas,"Liberation Mono",Menlo,Courier,monospace;color:#567482;background-color:#f3f6fa;border-radius:.3rem}.main-content pre{padding:.8rem;margin-top:0;margin-bottom:1rem;font:1rem Consolas,"Liberation Mono",Menlo,Courier,monospace;color:#567482;word-wrap:normal;background-color:#f3f6fa;border:solid 1px #dce6f0;border-radius:.3rem;line-height:1.45;overflow:auto}@media screen and (max-width:42em){.main-content pre{font-size:.9rem}}.main-content pre>code{padding:0;margin:0;color:#567482;word-break:normal;white-space:pre;background:0 0;border:0}@media screen and (max-width:42em){.main-content pre>code{font-size:.9rem}}.main-content pre code,.main-content pre tt{display:inline;max-width:initial;padding:0;margin:0;overflow:initial;line-height:inherit;word-wrap:normal;background-color:transparent;border:0}.main-content pre code:after,.main-content pre code:before,.main-content pre tt:after,.main-content pre tt:before{content:normal}.main-content ol,.main-content ul{margin-top:0}.main-content blockquote{padding:0 1rem;margin-left:0;color:#819198;border-left:.3rem solid #dce6f0;font-size:1.2rem}.main-content blockquote>:first-child{margin-top:0}.main-content blockquote>:last-child{margin-bottom:0}@media screen and (max-width:42em){.main-content blockquote{font-size:1.1rem}}.main-content table{width:100%;overflow:auto;word-break:normal;word-break:keep-all;-webkit-overflow-scrolling:touch;border-collapse:collapse;border-spacing:0;margin:1rem 0}.main-content table th{font-weight:700;background-color:#159957;color:#fff}.main-content table td,.main-content table th{padding:.5rem 1rem;border-bottom:1px solid #e9ebec;text-align:left}.main-content table tr:nth-child(odd){background-color:#f2f2f2}.main-content dl{padding:0}.main-content dl dt{padding:0;margin-top:1rem;font-size:1rem;font-weight:700}.main-content dl dd{padding:0;margin-bottom:1rem}.main-content hr{height:2px;padding:0;margin:1rem 0;background-color:#eff0f1;border:0}code span.kw { color: #a71d5d; font-weight: normal; } 
code span.dt { color: #795da3; } 
code span.dv { color: #0086b3; } 
code span.bn { color: #0086b3; } 
code span.fl { color: #0086b3; } 
code span.ch { color: #4070a0; } 
code span.st { color: #183691; } 
code span.co { color: #969896; font-style: italic; } 
code span.ot { color: #007020; } 
</style>





</head>

<body>




<section class="page-header">
<h1 class="title toc-ignore project-name">Bayesian Linear Regression</h1>
<h4 class="author project-author">庄亮亮</h4>
<h4 class="date project-date">2020/11/10</h4>
</section>



<section class="main-content">
<div class="sourceCode" id="cb1"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb1-1" title="1"><span class="co">#devtools::install_github(&quot;julianfaraway/brinla&quot;) #可能要翻墙</span></a>
<a class="sourceLine" id="cb1-2" title="2"><span class="kw">library</span>(INLA)</a>
<a class="sourceLine" id="cb1-3" title="3"><span class="kw">library</span>(brinla)</a>
<a class="sourceLine" id="cb1-4" title="4"><span class="kw">library</span>(ggplot2)</a>
<a class="sourceLine" id="cb1-5" title="5"><span class="kw">library</span>(GGally)</a>
<a class="sourceLine" id="cb1-6" title="6"><span class="kw">library</span>(tidyr)</a>
<a class="sourceLine" id="cb1-7" title="7"><span class="kw">library</span>(MASS)</a>
<a class="sourceLine" id="cb1-8" title="8"><span class="kw">library</span>(nlme)</a></code></pre></div>
<div id="bayesian-inference-for-linear-regression" class="section level1">
<h1>1. Bayesian inference for linear regression</h1>
<p><strong>空气污染数据</strong>：调查美国41个城市污染的决定因素。以<code>SO2</code>作为因变量，其他六个变量作为潜在解释变量。</p>
<p>在这些潜在解释变量中，有两个与人类生态相关<code>(pop,manuf)</code>，另外四个与气候相关<code>(negtemp,wind,precip,days)</code>。变量<code>negtemp</code>表示年平均气温的负值。在这里使用负值是因为所有的变量都是这样的，高的值表示一个不太事宜的环境。</p>
<div class="sourceCode" id="cb2"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb2-1" title="1"><span class="kw">data</span>(usair, <span class="dt">package =</span> <span class="st">&quot;brinla&quot;</span>)</a>
<a class="sourceLine" id="cb2-2" title="2">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">head</span>(usair), </a>
<a class="sourceLine" id="cb2-3" title="3">             <span class="dt">caption =</span> <span class="st">&#39;usair数据&#39;</span>, </a>
<a class="sourceLine" id="cb2-4" title="4">             <span class="dt">align=</span><span class="st">&#39;c&#39;</span>)</a></code></pre></div>
<table>
<caption>usair数据</caption>
<thead>
<tr class="header">
<th></th>
<th align="center">SO2</th>
<th align="center">negtemp</th>
<th align="center">manuf</th>
<th align="center">pop</th>
<th align="center">wind</th>
<th align="center">precip</th>
<th align="center">days</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>Phoenix</td>
<td align="center">10</td>
<td align="center">-70.3</td>
<td align="center">213</td>
<td align="center">582</td>
<td align="center">6.0</td>
<td align="center">7.05</td>
<td align="center">36</td>
</tr>
<tr class="even">
<td>Little Rock</td>
<td align="center">13</td>
<td align="center">-61.0</td>
<td align="center">91</td>
<td align="center">132</td>
<td align="center">8.2</td>
<td align="center">48.52</td>
<td align="center">100</td>
</tr>
<tr class="odd">
<td>San Francisco</td>
<td align="center">12</td>
<td align="center">-56.7</td>
<td align="center">453</td>
<td align="center">716</td>
<td align="center">8.7</td>
<td align="center">20.66</td>
<td align="center">67</td>
</tr>
<tr class="even">
<td>Denver</td>
<td align="center">17</td>
<td align="center">-51.9</td>
<td align="center">454</td>
<td align="center">515</td>
<td align="center">9.0</td>
<td align="center">12.95</td>
<td align="center">86</td>
</tr>
<tr class="odd">
<td>Hartford</td>
<td align="center">56</td>
<td align="center">-49.1</td>
<td align="center">412</td>
<td align="center">158</td>
<td align="center">9.0</td>
<td align="center">43.37</td>
<td align="center">127</td>
</tr>
<tr class="even">
<td>Wilmington</td>
<td align="center">36</td>
<td align="center">-54.0</td>
<td align="center">80</td>
<td align="center">80</td>
<td align="center">9.0</td>
<td align="center">40.25</td>
<td align="center">114</td>
</tr>
</tbody>
</table>
<div class="sourceCode" id="cb3"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb3-1" title="1">pairs.chart &lt;-<span class="st"> </span><span class="kw">ggpairs</span>(usair[,<span class="op">-</span><span class="dv">1</span>], <span class="dt">lower =</span> <span class="kw">list</span>(<span class="dt">continuous =</span> <span class="st">&quot;cor&quot;</span>), <span class="dt">upper =</span> <span class="kw">list</span>(<span class="dt">continuous =</span> <span class="st">&quot;points&quot;</span>, <span class="dt">combo =</span> <span class="st">&quot;dot&quot;</span>)) <span class="op">+</span><span class="st"> </span>ggplot2<span class="op">::</span><span class="kw">theme</span>(<span class="dt">axis.text =</span> <span class="kw">element_text</span>(<span class="dt">size =</span> <span class="dv">6</span>)) </a>
<a class="sourceLine" id="cb3-2" title="2">pairs.chart  </a></code></pre></div>
<p><img src="" /><!-- --></p>
<p><code>Manuf</code>与<code>Pop</code>高度相关:<code>r=0.955</code>；我们在模型中只需保留一个。</p>
<div id="frequentist-approach" class="section level2">
<h2>1.1 Frequentist approach</h2>
<div class="sourceCode" id="cb4"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb4-1" title="1">usair.formula1 &lt;-<span class="st">  </span>SO2 <span class="op">~</span><span class="st"> </span>negtemp <span class="op">+</span><span class="st"> </span>manuf <span class="op">+</span><span class="st"> </span>wind <span class="op">+</span><span class="st"> </span>precip <span class="op">+</span><span class="st"> </span>days</a>
<a class="sourceLine" id="cb4-2" title="2">usair.lm1 &lt;-<span class="st"> </span><span class="kw">lm</span>(usair.formula1, <span class="dt">data =</span> usair)</a>
<a class="sourceLine" id="cb4-3" title="3">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">round</span>(<span class="kw">coef</span>(<span class="kw">summary</span>(usair.lm1)), <span class="dv">4</span>), </a>
<a class="sourceLine" id="cb4-4" title="4">             <span class="dt">caption =</span> <span class="st">&#39;拟合情况汇总&#39;</span>, </a>
<a class="sourceLine" id="cb4-5" title="5">             <span class="dt">align=</span><span class="st">&#39;c&#39;</span>)</a></code></pre></div>
<table>
<caption>拟合情况汇总</caption>
<thead>
<tr class="header">
<th></th>
<th align="center">Estimate</th>
<th align="center">Std. Error</th>
<th align="center">t value</th>
<th align="center">Pr(&gt;|t|)</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>(Intercept)</td>
<td align="center">135.7714</td>
<td align="center">50.0610</td>
<td align="center">2.7121</td>
<td align="center">0.0103</td>
</tr>
<tr class="even">
<td>negtemp</td>
<td align="center">1.7714</td>
<td align="center">0.6366</td>
<td align="center">2.7824</td>
<td align="center">0.0086</td>
</tr>
<tr class="odd">
<td>manuf</td>
<td align="center">0.0256</td>
<td align="center">0.0046</td>
<td align="center">5.5544</td>
<td align="center">0.0000</td>
</tr>
<tr class="even">
<td>wind</td>
<td align="center">-3.7379</td>
<td align="center">1.9444</td>
<td align="center">-1.9224</td>
<td align="center">0.0627</td>
</tr>
<tr class="odd">
<td>precip</td>
<td align="center">0.6259</td>
<td align="center">0.3885</td>
<td align="center">1.6111</td>
<td align="center">0.1161</td>
</tr>
<tr class="even">
<td>days</td>
<td align="center">-0.0571</td>
<td align="center">0.1748</td>
<td align="center">-0.3265</td>
<td align="center">0.7460</td>
</tr>
</tbody>
</table>
<p>结果表明，<code>negtemp</code>和<code>manuf</code>是重要的解释变量，而<code>wind</code>,<code>precip</code>和<code>days</code>不是。</p>
<ul>
<li>标准残差</li>
</ul>
<div class="sourceCode" id="cb5"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb5-1" title="1"><span class="kw">round</span>(<span class="kw">summary</span>(usair.lm1)<span class="op">$</span>sigma, <span class="dv">4</span>)</a></code></pre></div>
<pre><code>## [1] 15.79</code></pre>
</div>
<div id="bayesian-regression-with-inla" class="section level2">
<h2>1.2 Bayesian Regression with INLA</h2>
<p>默认情况下： <span class="math display">\[\begin{equation*}
\beta_{j} \sim N\left(0,10^{6}\right), \quad j=0, \ldots, p
\end{equation*}\]</span></p>
<p><span class="math display">\[\begin{equation*}
\log (\tau) \sim \log \operatorname{Gamma}\left(1,10^{-5}\right)
\end{equation*}\]</span></p>
<div class="sourceCode" id="cb7"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb7-1" title="1">usair.inla1 &lt;-<span class="st"> </span><span class="kw">inla</span>(usair.formula1, <span class="dt">data =</span> usair, <span class="dt">control.compute =</span> <span class="kw">list</span>(<span class="dt">dic =</span> <span class="ot">TRUE</span>, <span class="dt">cpo =</span> <span class="ot">TRUE</span>))</a></code></pre></div>
<div class="sourceCode" id="cb8"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb8-1" title="1">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">round</span>(usair.inla1<span class="op">$</span>summary.fixed, <span class="dv">4</span>), </a>
<a class="sourceLine" id="cb8-2" title="2">             <span class="dt">caption =</span> <span class="st">&#39;固定效应信息&#39;</span>, </a>
<a class="sourceLine" id="cb8-3" title="3">             <span class="dt">align=</span><span class="st">&#39;c&#39;</span>)</a></code></pre></div>
<table>
<caption>固定效应信息</caption>
<thead>
<tr class="header">
<th></th>
<th align="center">mean</th>
<th align="center">sd</th>
<th align="center">0.025quant</th>
<th align="center">0.5quant</th>
<th align="center">0.975quant</th>
<th align="center">mode</th>
<th align="center">kld</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>(Intercept)</td>
<td align="center">135.4904</td>
<td align="center">49.8882</td>
<td align="center">36.9784</td>
<td align="center">135.4963</td>
<td align="center">233.8821</td>
<td align="center">135.5120</td>
<td align="center">0</td>
</tr>
<tr class="even">
<td>negtemp</td>
<td align="center">1.7690</td>
<td align="center">0.6347</td>
<td align="center">0.5157</td>
<td align="center">1.7690</td>
<td align="center">3.0209</td>
<td align="center">1.7692</td>
<td align="center">0</td>
</tr>
<tr class="odd">
<td>manuf</td>
<td align="center">0.0256</td>
<td align="center">0.0046</td>
<td align="center">0.0165</td>
<td align="center">0.0256</td>
<td align="center">0.0346</td>
<td align="center">0.0256</td>
<td align="center">0</td>
</tr>
<tr class="even">
<td>wind</td>
<td align="center">-3.7230</td>
<td align="center">1.9357</td>
<td align="center">-7.5432</td>
<td align="center">-3.7234</td>
<td align="center">0.0964</td>
<td align="center">-3.7241</td>
<td align="center">0</td>
</tr>
<tr class="odd">
<td>precip</td>
<td align="center">0.6249</td>
<td align="center">0.3874</td>
<td align="center">-0.1401</td>
<td align="center">0.6249</td>
<td align="center">1.3890</td>
<td align="center">0.6249</td>
<td align="center">0</td>
</tr>
<tr class="even">
<td>days</td>
<td align="center">-0.0567</td>
<td align="center">0.1743</td>
<td align="center">-0.4007</td>
<td align="center">-0.0567</td>
<td align="center">0.2872</td>
<td align="center">-0.0567</td>
<td align="center">0</td>
</tr>
</tbody>
</table>
<div class="sourceCode" id="cb9"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb9-1" title="1">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">round</span>(usair.inla1<span class="op">$</span>summary.hyperpar, <span class="dv">4</span>), </a>
<a class="sourceLine" id="cb9-2" title="2">             <span class="dt">caption =</span> <span class="st">&#39;超参数信息&#39;</span>, </a>
<a class="sourceLine" id="cb9-3" title="3">             <span class="dt">align=</span><span class="st">&#39;c&#39;</span>)</a></code></pre></div>
<table>
<caption>超参数信息</caption>
<thead>
<tr class="header">
<th></th>
<th align="center">mean</th>
<th align="center">sd</th>
<th align="center">0.025quant</th>
<th align="center">0.5quant</th>
<th align="center">0.975quant</th>
<th align="center">mode</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>Precision for the Gaussian observations</td>
<td align="center">0.0042</td>
<td align="center">0.001</td>
<td align="center">0.0025</td>
<td align="center">0.0042</td>
<td align="center">0.0064</td>
<td align="center">0.004</td>
</tr>
</tbody>
</table>
<div class="sourceCode" id="cb10"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb10-1" title="1"><span class="kw">summary</span>(usair.inla1)</a></code></pre></div>
<pre><code>## 
## Call:
##    c(&quot;inla(formula = usair.formula1, data = usair, control.compute = 
##    list(dic = TRUE, &quot;, &quot; cpo = TRUE))&quot;) 
## Time used:
##     Pre = 1.26, Running = 0.386, Post = 0.125, Total = 1.77 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) 135.490 49.888     36.978  135.496    233.882 135.512   0
## negtemp       1.769  0.635      0.516    1.769      3.021   1.769   0
## manuf         0.026  0.005      0.017    0.026      0.035   0.026   0
## wind         -3.723  1.936     -7.543   -3.723      0.096  -3.724   0
## precip        0.625  0.387     -0.140    0.625      1.389   0.625   0
## days         -0.057  0.174     -0.401   -0.057      0.287  -0.057   0
## 
## Model hyperparameters:
##                                          mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.004 0.001      0.003    0.004
##                                         0.975quant  mode
## Precision for the Gaussian observations      0.006 0.004
## 
## Expected number of effective parameters(stdev): 6.00(0.001)
## Number of equivalent replicates : 6.84 
## 
## Deviance Information Criterion (DIC) ...............: 350.76
## Deviance Information Criterion (DIC, saturated) ....: 51.32
## Effective number of parameters .....................: 7.21
## 
## Marginal log-Likelihood:  -208.73 
## CPO and PIT are computed
## 
## Posterior marginals for the linear predictor and
##  the fitted values are computed</code></pre>
<ul>
<li>计算<span class="math inline">\(\sigma\)</span>的后验估计值</li>
</ul>
<p>默认情况下，inla对象输出的是后验信息的精确参数<span class="math inline">\(\tau\)</span>。然而常常我们对后验均值<span class="math inline">\(\sigma\)</span>感兴趣。</p>
<p><code>bri.hyperpar.summary</code>用于生成超参数<span class="math inline">\(\sigma\)</span>的汇总统计信息。</p>
<div class="sourceCode" id="cb12"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb12-1" title="1">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">round</span>(<span class="kw">bri.hyperpar.summary</span>(usair.inla1), <span class="dv">4</span>), </a>
<a class="sourceLine" id="cb12-2" title="2">             <span class="dt">caption =</span> <span class="st">&#39;超参数信息&#39;</span>, </a>
<a class="sourceLine" id="cb12-3" title="3">             <span class="dt">align=</span><span class="st">&#39;c&#39;</span>)</a></code></pre></div>
<table>
<caption>超参数信息</caption>
<thead>
<tr class="header">
<th></th>
<th align="center">mean</th>
<th align="center">sd</th>
<th align="center">q0.025</th>
<th align="center">q0.5</th>
<th align="center">q0.975</th>
<th align="center">mode</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>SD for the Gaussian observations</td>
<td align="center">15.67</td>
<td align="center">1.865</td>
<td align="center">12.53</td>
<td align="center">15.49</td>
<td align="center">19.84</td>
<td align="center">15.14</td>
</tr>
</tbody>
</table>
<ul>
<li>绘制<span class="math inline">\(\sigma\)</span>的后验密度</li>
</ul>
<p>使用函数<code>bri.hyperpar.plot</code>。</p>
<div class="sourceCode" id="cb13"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb13-1" title="1"><span class="kw">bri.hyperpar.plot</span>(usair.inla1)</a></code></pre></div>
<p><img src="" /><!-- --></p>
<p><span class="math inline">\(\sigma\)</span>有略微右偏的后验分布。</p>
<ul>
<li>改变先验</li>
</ul>
<p>假设 <span class="math inline">\(\beta_0 \sim N(100,100)\)</span>, <span class="math inline">\(\beta_{negtemp} \sim N(2,1)\)</span> 和 <span class="math inline">\(\beta_{wind} \sim N(−3,1)\)</span>。</p>
<p>如果想假设<span class="math inline">\(\tau\)</span>服从对数正态分布(对数<span class="math inline">\(\tau\)</span>服从正态分布)，可以使用选项<code>control.family</code>来指定。</p>
<div class="sourceCode" id="cb14"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb14-1" title="1">usair.inla2 &lt;-<span class="st"> </span><span class="kw">inla</span>(usair.formula1, <span class="dt">data =</span> usair, <span class="dt">control.compute =</span> <span class="kw">list</span>(<span class="dt">dic =</span> <span class="ot">TRUE</span>, <span class="dt">cpo =</span><span class="ot">TRUE</span>),</a>
<a class="sourceLine" id="cb14-2" title="2">    <span class="dt">control.fixed =</span> <span class="kw">list</span>(<span class="dt">mean.intercept =</span> <span class="dv">100</span>, <span class="dt">prec.intercept =</span> <span class="dv">10</span><span class="op">^</span>(<span class="op">-</span><span class="dv">2</span>), <span class="dt">mean =</span> <span class="kw">list</span>(<span class="dt">negtemp =</span> <span class="dv">2</span>, <span class="dt">wind =</span> <span class="dv">-3</span>, <span class="dt">default =</span><span class="dv">0</span>), <span class="dt">prec =</span> <span class="dv">1</span>), </a>
<a class="sourceLine" id="cb14-3" title="3">    <span class="dt">control.family =</span> <span class="kw">list</span>(<span class="dt">hyper =</span> <span class="kw">list</span>(<span class="dt">prec =</span> <span class="kw">list</span>(<span class="dt">prior=</span><span class="st">&quot;gaussian&quot;</span>, <span class="dt">param =</span><span class="kw">c</span>(<span class="dv">0</span>,<span class="dv">1</span>)))))</a>
<a class="sourceLine" id="cb14-4" title="4"><span class="kw">summary</span>(usair.inla2)</a></code></pre></div>
<pre><code>## 
## Call:
##    c(&quot;inla(formula = usair.formula1, data = usair, control.compute = 
##    list(dic = TRUE, &quot;, &quot; cpo = TRUE), control.family = list(hyper = 
##    list(prec = list(prior = \&quot;gaussian\&quot;, &quot;, &quot; param = c(0, 1)))), 
##    control.fixed = list(mean.intercept = 100, &quot;, &quot; prec.intercept = 
##    10^(-2), mean = list(negtemp = 2, wind = -3, &quot;, &quot; default = 0), prec = 
##    1))&quot;) 
## Time used:
##     Pre = 0.586, Running = 0.4, Post = 0.119, Total = 1.1 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) 102.395 9.562     83.615  102.397    121.149 102.401   0
## negtemp       1.384 0.212      0.966    1.385      1.801   1.385   0
## manuf         0.025 0.004      0.018    0.025      0.033   0.025   0
## wind         -2.952 0.803     -4.529   -2.952     -1.376  -2.952   0
## precip        0.441 0.246     -0.045    0.441      0.925   0.442   0
## days          0.041 0.090     -0.136    0.040      0.218   0.040   0
## 
## Model hyperparameters:
##                                          mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.005 0.001      0.003    0.005
##                                         0.975quant  mode
## Precision for the Gaussian observations      0.008 0.005
## 
## Expected number of effective parameters(stdev): 4.33(0.074)
## Number of equivalent replicates : 9.48 
## 
## Deviance Information Criterion (DIC) ...............: 347.57
## Deviance Information Criterion (DIC, saturated) ....: 57.38
## Effective number of parameters .....................: 5.26
## 
## Marginal log-Likelihood:  -197.50 
## CPO and PIT are computed
## 
## Posterior marginals for the linear predictor and
##  the fitted values are computed</code></pre>
<hr />
</div>
</div>
<div id="prediction" class="section level1">
<h1>2. Prediction</h1>
<div id="the-prediction-in-lm" class="section level2">
<h2>2.1 The prediction in lm</h2>
<div class="sourceCode" id="cb16"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb16-1" title="1"><span class="co"># new observations</span></a>
<a class="sourceLine" id="cb16-2" title="2">new.data &lt;-<span class="st"> </span><span class="kw">data.frame</span>(</a>
<a class="sourceLine" id="cb16-3" title="3">  <span class="dt">negtemp =</span> <span class="kw">c</span>(<span class="op">-</span><span class="dv">50</span>, <span class="dv">-60</span>, <span class="dv">-40</span>),</a>
<a class="sourceLine" id="cb16-4" title="4">  <span class="dt">manuf =</span> <span class="kw">c</span>(<span class="dv">150</span>, <span class="dv">100</span>, <span class="dv">400</span>), </a>
<a class="sourceLine" id="cb16-5" title="5">  <span class="dt">pop =</span> <span class="kw">c</span>(<span class="dv">200</span>, <span class="dv">100</span>, <span class="dv">300</span>), </a>
<a class="sourceLine" id="cb16-6" title="6">  <span class="dt">wind =</span> <span class="kw">c</span>(<span class="dv">6</span>, <span class="dv">7</span>, <span class="dv">8</span>), </a>
<a class="sourceLine" id="cb16-7" title="7">  <span class="dt">precip =</span> <span class="kw">c</span>(<span class="dv">10</span>, <span class="dv">30</span>, <span class="dv">20</span>), </a>
<a class="sourceLine" id="cb16-8" title="8">  <span class="dt">days =</span> <span class="kw">c</span>(<span class="dv">20</span>, <span class="dv">100</span>, <span class="dv">40</span>))</a>
<a class="sourceLine" id="cb16-9" title="9">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">head</span>(new.data))</a></code></pre></div>
<table>
<thead>
<tr class="header">
<th align="right">negtemp</th>
<th align="right">manuf</th>
<th align="right">pop</th>
<th align="right">wind</th>
<th align="right">precip</th>
<th align="right">days</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td align="right">-50</td>
<td align="right">150</td>
<td align="right">200</td>
<td align="right">6</td>
<td align="right">10</td>
<td align="right">20</td>
</tr>
<tr class="even">
<td align="right">-60</td>
<td align="right">100</td>
<td align="right">100</td>
<td align="right">7</td>
<td align="right">30</td>
<td align="right">100</td>
</tr>
<tr class="odd">
<td align="right">-40</td>
<td align="right">400</td>
<td align="right">300</td>
<td align="right">8</td>
<td align="right">20</td>
<td align="right">40</td>
</tr>
</tbody>
</table>
<div class="sourceCode" id="cb17"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb17-1" title="1"><span class="kw">predict</span>(usair.lm1, new.data, <span class="dt">se.fit =</span> <span class="ot">TRUE</span>)</a></code></pre></div>
<pre><code>## $fit
##     1     2     3 
## 33.73 18.95 55.48 
## 
## $se.fit
##      1      2      3 
## 14.937  5.329 17.639 
## 
## $df
## [1] 35
## 
## $residual.scale
## [1] 15.79</code></pre>
<p>输出：<code>预测向量(\$fit)</code>、<code>预测均值的标准误差向量(\$se.fit)</code>、<code>残差的自由度(\$df)</code>和<code>残差标准差($residual.scale)</code>。</p>
</div>
<div id="the-prediction-in-inla" class="section level2">
<h2>2.2 The prediction in INLA</h2>
<p>在INLA中，与lm的函数预测不同。预测可以作为模型拟合的一部分。由于预测和拟合一个有缺失数据的模型是一样的，我们需要设置响应变量<code>“y[i]=NA”</code>表示我们想要预测的值。</p>
<div class="sourceCode" id="cb19"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb19-1" title="1">usair.combined &lt;-<span class="st"> </span><span class="kw">rbind</span>(usair, <span class="kw">data.frame</span>(<span class="dt">SO2 =</span> <span class="kw">c</span>(<span class="ot">NA</span>, <span class="ot">NA</span>, <span class="ot">NA</span>), new.data))</a>
<a class="sourceLine" id="cb19-2" title="2">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">tail</span>(usair.combined))</a></code></pre></div>
<table>
<thead>
<tr class="header">
<th></th>
<th align="right">SO2</th>
<th align="right">negtemp</th>
<th align="right">manuf</th>
<th align="right">pop</th>
<th align="right">wind</th>
<th align="right">precip</th>
<th align="right">days</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>Seattle</td>
<td align="right">29</td>
<td align="right">-51.1</td>
<td align="right">379</td>
<td align="right">531</td>
<td align="right">9.4</td>
<td align="right">38.79</td>
<td align="right">164</td>
</tr>
<tr class="even">
<td>Charleston</td>
<td align="right">31</td>
<td align="right">-55.2</td>
<td align="right">35</td>
<td align="right">71</td>
<td align="right">6.5</td>
<td align="right">40.75</td>
<td align="right">148</td>
</tr>
<tr class="odd">
<td>Milwaukee</td>
<td align="right">16</td>
<td align="right">-45.7</td>
<td align="right">569</td>
<td align="right">717</td>
<td align="right">11.8</td>
<td align="right">29.07</td>
<td align="right">123</td>
</tr>
<tr class="even">
<td>1</td>
<td align="right">NA</td>
<td align="right">-50.0</td>
<td align="right">150</td>
<td align="right">200</td>
<td align="right">6.0</td>
<td align="right">10.00</td>
<td align="right">20</td>
</tr>
<tr class="odd">
<td>2</td>
<td align="right">NA</td>
<td align="right">-60.0</td>
<td align="right">100</td>
<td align="right">100</td>
<td align="right">7.0</td>
<td align="right">30.00</td>
<td align="right">100</td>
</tr>
<tr class="even">
<td>3</td>
<td align="right">NA</td>
<td align="right">-40.0</td>
<td align="right">400</td>
<td align="right">300</td>
<td align="right">8.0</td>
<td align="right">20.00</td>
<td align="right">40</td>
</tr>
</tbody>
</table>
<div class="sourceCode" id="cb20"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb20-1" title="1">usair.link &lt;-<span class="st"> </span><span class="kw">c</span>(<span class="kw">rep</span>(<span class="ot">NA</span>, <span class="kw">nrow</span>(usair)), <span class="kw">rep</span>(<span class="dv">1</span>, <span class="kw">nrow</span>(new.data)))</a>
<a class="sourceLine" id="cb20-2" title="2"><span class="kw">tail</span>(usair.link)</a></code></pre></div>
<pre><code>## [1] NA NA NA  1  1  1</code></pre>
<div class="sourceCode" id="cb22"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb22-1" title="1">usair.inla1.pred &lt;-<span class="st"> </span><span class="kw">inla</span>(usair.formula1, <span class="dt">data =</span> usair.combined, <span class="dt">control.predictor =</span> <span class="kw">list</span>(<span class="dt">link =</span> usair.link))</a>
<a class="sourceLine" id="cb22-2" title="2"></a>
<a class="sourceLine" id="cb22-3" title="3">knitr<span class="op">::</span><span class="kw">kable</span>(usair.inla1.pred<span class="op">$</span>summary.fitted.values[(<span class="kw">nrow</span>(usair)<span class="op">+</span><span class="dv">1</span>)<span class="op">:</span><span class="kw">nrow</span>(usair.combined),], <span class="dt">caption =</span> <span class="st">&#39;拟合情况&#39;</span>, <span class="dt">align=</span><span class="st">&#39;c&#39;</span>)</a></code></pre></div>
<table>
<caption>拟合情况</caption>
<thead>
<tr class="header">
<th></th>
<th align="center">mean</th>
<th align="center">sd</th>
<th align="center">0.025quant</th>
<th align="center">0.5quant</th>
<th align="center">0.975quant</th>
<th align="center">mode</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>fitted.Predictor.42</td>
<td align="center">33.65</td>
<td align="center">14.88</td>
<td align="center">4.316</td>
<td align="center">33.66</td>
<td align="center">63.00</td>
<td align="center">33.66</td>
</tr>
<tr class="even">
<td>fitted.Predictor.43</td>
<td align="center">18.93</td>
<td align="center">5.31</td>
<td align="center">8.461</td>
<td align="center">18.93</td>
<td align="center">29.40</td>
<td align="center">18.93</td>
</tr>
<tr class="odd">
<td>fitted.Predictor.44</td>
<td align="center">55.41</td>
<td align="center">17.58</td>
<td align="center">20.753</td>
<td align="center">55.41</td>
<td align="center">90.07</td>
<td align="center">55.41</td>
</tr>
</tbody>
</table>
<hr />
</div>
</div>
<div id="model-selection-and-checking" class="section level1">
<h1>3. Model selection and checking</h1>
<div id="dic" class="section level2">
<h2>3.1 DIC</h2>
<ul>
<li>AIC</li>
</ul>
<div class="sourceCode" id="cb23"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb23-1" title="1">usair.step &lt;-<span class="st"> </span><span class="kw">stepAIC</span>(usair.lm1, <span class="dt">trace =</span> <span class="ot">FALSE</span>)</a>
<a class="sourceLine" id="cb23-2" title="2">usair.step<span class="op">$</span>anova</a></code></pre></div>
<pre><code>## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## SO2 ~ negtemp + manuf + wind + precip + days
## 
## Final Model:
## SO2 ~ negtemp + manuf + wind + precip
## 
## 
##     Step Df Deviance Resid. Df Resid. Dev   AIC
## 1                           35       8726 231.8
## 2 - days  1    26.57        36       8753 229.9</code></pre>
<div class="sourceCode" id="cb25"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb25-1" title="1"><span class="co"># Final multiple regression model</span></a>
<a class="sourceLine" id="cb25-2" title="2">usair.formula2 &lt;-<span class="st"> </span>SO2 <span class="op">~</span><span class="st"> </span>negtemp <span class="op">+</span><span class="st"> </span>manuf <span class="op">+</span><span class="st"> </span>wind <span class="op">+</span><span class="st"> </span>precip</a>
<a class="sourceLine" id="cb25-3" title="3">usair.lm2 &lt;-<span class="st"> </span><span class="kw">lm</span>(usair.formula2, <span class="dt">data =</span> usair)</a>
<a class="sourceLine" id="cb25-4" title="4">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">round</span>(<span class="kw">coef</span>(<span class="kw">summary</span>(usair.lm2)), <span class="dv">4</span>))</a></code></pre></div>
<table>
<thead>
<tr class="header">
<th></th>
<th align="right">Estimate</th>
<th align="right">Std. Error</th>
<th align="right">t value</th>
<th align="right">Pr(&gt;|t|)</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>(Intercept)</td>
<td align="right">123.1183</td>
<td align="right">31.2907</td>
<td align="right">3.935</td>
<td align="right">0.0004</td>
</tr>
<tr class="even">
<td>negtemp</td>
<td align="right">1.6114</td>
<td align="right">0.4014</td>
<td align="right">4.015</td>
<td align="right">0.0003</td>
</tr>
<tr class="odd">
<td>manuf</td>
<td align="right">0.0255</td>
<td align="right">0.0045</td>
<td align="right">5.615</td>
<td align="right">0.0000</td>
</tr>
<tr class="even">
<td>wind</td>
<td align="right">-3.6302</td>
<td align="right">1.8923</td>
<td align="right">-1.918</td>
<td align="right">0.0630</td>
</tr>
<tr class="odd">
<td>precip</td>
<td align="right">0.5242</td>
<td align="right">0.2294</td>
<td align="right">2.285</td>
<td align="right">0.0283</td>
</tr>
</tbody>
</table>
<ul>
<li>DIC</li>
</ul>
<p>在贝叶斯分析中，DIC是AIC的推广，是最常用的贝叶斯模型比较方法之一，定义为拟合优度度量和模型复杂度度量的总和</p>
<div class="sourceCode" id="cb26"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb26-1" title="1">usair.inla3 &lt;-<span class="st"> </span><span class="kw">inla</span>(usair.formula2, <span class="dt">data =</span> usair, <span class="dt">control.compute =</span> <span class="kw">list</span>(<span class="dt">dic =</span> <span class="ot">TRUE</span>, <span class="dt">cpo =</span> <span class="ot">TRUE</span>))</a>
<a class="sourceLine" id="cb26-2" title="2"><span class="kw">c</span>(usair.inla1<span class="op">$</span>dic<span class="op">$</span>dic, usair.inla3<span class="op">$</span>dic<span class="op">$</span>dic)</a></code></pre></div>
<pre><code>## [1] 350.8 348.7</code></pre>
<p>最佳模型：所有子集回归中最小的DIC。</p>
</div>
<div id="posterior-predictive-checking" class="section level2">
<h2>3.2 Posterior Predictive Checking</h2>
<p>在贝叶斯分析中，模型评估通常是</p>
<ol style="list-style-type: decimal">
<li><p>基于后验预测检查；</p></li>
<li><p>留一交叉验证预测检查。</p></li>
</ol>
<div id="基于后验预测检查" class="section level3">
<h3>3.2.1 基于后验预测检查</h3>
<p>后验预测p值可以通过R函数<code>INLA .pmarginal</code>得到</p>
<div class="sourceCode" id="cb28"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb28-1" title="1">usair.inla3.pred &lt;-<span class="st"> </span><span class="kw">inla</span>(usair.formula2, <span class="dt">data =</span> usair, <span class="dt">control.predictor =</span> <span class="kw">list</span>(<span class="dt">link =</span> <span class="dv">1</span>, <span class="dt">compute =</span> <span class="ot">TRUE</span>))</a>
<a class="sourceLine" id="cb28-2" title="2">post.predicted.pval &lt;-<span class="st"> </span><span class="kw">vector</span>(<span class="dt">mode =</span> <span class="st">&quot;numeric&quot;</span>, <span class="dt">length =</span> <span class="kw">nrow</span>(usair))</a>
<a class="sourceLine" id="cb28-3" title="3"><span class="cf">for</span>(i <span class="cf">in</span> (<span class="dv">1</span><span class="op">:</span><span class="kw">nrow</span>(usair))) {</a>
<a class="sourceLine" id="cb28-4" title="4">  post.predicted.pval[i] &lt;-<span class="st"> </span><span class="kw">inla.pmarginal</span>(<span class="dt">q=</span>usair<span class="op">$</span>SO2[i], <span class="dt">marginal =</span> usair.inla3.pred<span class="op">$</span>marginals.fitted.values[[i]])</a>
<a class="sourceLine" id="cb28-5" title="5">}</a>
<a class="sourceLine" id="cb28-6" title="6"><span class="kw">hist</span>(post.predicted.pval, <span class="dt">main=</span><span class="st">&quot;&quot;</span>, <span class="dt">breaks =</span> <span class="dv">10</span>, <span class="dt">xlab=</span><span class="st">&quot;Posterior predictive p-value&quot;</span>)</a></code></pre></div>
<p><img src="" /><!-- --></p>
<p>很多后验预测p值都接近于0或1。然而，解释后验预测p值的一个<strong>缺点</strong>是，即使数据来自真实的模型，它们也不能具有均匀分布。 因此，后验预测p值的图并不令人满意，我们希望使用其他模型评估方法进一步检验该模型。</p>
</div>
<div id="留一交叉验证预测检查" class="section level3">
<h3>3.2.2 留一交叉验证预测检查</h3>
<p>评价模型的优劣的两个量：</p>
<ol style="list-style-type: decimal">
<li>conditional predictive ordinate (CPO)</li>
</ol>
<p><span class="math display">\[
    CPO_{i} =p\left(y_{i} \mid y_{-i}\right) 
\]</span></p>
<ol start="2" style="list-style-type: decimal">
<li>probability integral transform (PIT)</li>
</ol>
<p><span class="math display">\[
    \mathrm{PIT}_{i} =p\left(y_{i}^{*} \leq y_{i} \mid \mathbf{y}_{-i}\right)
\]</span></p>
<p>在INLA中有针对潜在问题的内部检查，这些问题出现在<code>usair.inla3$cpo$failure</code>中。它是一个向量，每个观测值都包含0或1。当值为1时，表示CPO或PIT的估计对于相应的观测是不可靠的。在我们的例子中，我们可以通过以下方法检查是否存在故障：</p>
<div class="sourceCode" id="cb29"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb29-1" title="1"><span class="kw">sum</span>(usair.inla3<span class="op">$</span>cpo<span class="op">$</span>failure)</a></code></pre></div>
<pre><code>## [1] 0</code></pre>
<p>因此，在<code>usair.inla3</code>中不存在CPOs和PITs的计算问题。</p>
<div class="sourceCode" id="cb31"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb31-1" title="1"><span class="kw">par</span>(<span class="dt">mfrow =</span> <span class="kw">c</span>(<span class="dv">1</span>, <span class="dv">2</span>))</a>
<a class="sourceLine" id="cb31-2" title="2"><span class="kw">hist</span>(usair.inla3<span class="op">$</span>cpo<span class="op">$</span>pit, <span class="dt">main=</span><span class="st">&quot;&quot;</span>, <span class="dt">breaks =</span> <span class="dv">10</span>, <span class="dt">xlab =</span> <span class="st">&quot;PIT&quot;</span>)</a>
<a class="sourceLine" id="cb31-3" title="3"><span class="kw">qqplot</span>(<span class="kw">qunif</span>(<span class="kw">ppoints</span>(<span class="kw">length</span>(usair.inla3<span class="op">$</span>cpo<span class="op">$</span>pit))), usair.inla3<span class="op">$</span>cpo<span class="op">$</span>pit, <span class="dt">main =</span> <span class="st">&quot;Q-Q plot for Unif(0,1)&quot;</span>, <span class="dt">xlab =</span> <span class="st">&quot;Theoretical Quantiles&quot;</span>, <span class="dt">ylab =</span> <span class="st">&quot;Sample Quantiles&quot;</span>)</a>
<a class="sourceLine" id="cb31-4" title="4"><span class="kw">qqline</span>(usair.inla3<span class="op">$</span>cpo<span class="op">$</span>pit, <span class="dt">distribution =</span> <span class="cf">function</span>(p) <span class="kw">qunif</span>(p), <span class="dt">prob =</span> <span class="kw">c</span>(<span class="fl">0.1</span>, <span class="fl">0.9</span>))</a></code></pre></div>
<p><img src="" /><!-- --></p>
<p>PITs的分布接近均匀分布，表明模型对数据的拟合较为合理。值得注意的是，PITs直方图比对应的后验预测直方图更接近均匀分布。</p>
<ul>
<li>Compare LPML for the full model and reduce model</li>
</ul>
<p>如果我们把所有<span class="math inline">\(CPO\)</span>值的乘积看作一个“伪边际似然”，这就给出了一个交叉验证的拟合度度量。Geisser和Eddy(1979)提出的对数拟边际似然(<span class="math inline">\(LPML\)</span>)： <span class="math display">\[
L P M L=\log \left\{\prod_{i=1}^{n} p\left(y_{i} \mid \mathbf{y}_{-i}\right)\right\}=\sum_{i=1}^{n} \log p\left(y_{i} \mid \mathbf{y}_{-i}\right)=\sum_{i=1}^{n} \log \mathrm{CPO}_{i}
\]</span></p>
<div class="sourceCode" id="cb32"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb32-1" title="1">LPML1 &lt;-<span class="st"> </span><span class="kw">sum</span>(<span class="kw">log</span>(usair.inla1<span class="op">$</span>cpo<span class="op">$</span>cpo))</a>
<a class="sourceLine" id="cb32-2" title="2">LPML3 &lt;-<span class="st"> </span><span class="kw">sum</span>(<span class="kw">log</span>(usair.inla3<span class="op">$</span>cpo<span class="op">$</span>cpo))</a>
<a class="sourceLine" id="cb32-3" title="3"><span class="kw">c</span>(LPML1, LPML3)</a></code></pre></div>
<pre><code>## [1] -177.4 -176.1</code></pre>
<p>简化模型的LPML比完整模型的LPML大，这表明简化模型是首选的。</p>
<ul>
<li>CPOs的逐点比较</li>
</ul>
<div class="sourceCode" id="cb34"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb34-1" title="1"><span class="kw">par</span>(<span class="dt">mfrow=</span><span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">1</span>))</a>
<a class="sourceLine" id="cb34-2" title="2"><span class="kw">plot</span>(usair.inla1<span class="op">$</span>cpo<span class="op">$</span>cpo, usair.inla3<span class="op">$</span>cpo<span class="op">$</span>cpo, <span class="dt">xlab=</span><span class="st">&quot;CPO for the full model&quot;</span>, <span class="dt">ylab=</span><span class="st">&quot;CPO for the reduced model&quot;</span>)</a>
<a class="sourceLine" id="cb34-3" title="3"><span class="kw">abline</span>(<span class="dv">0</span>,<span class="dv">1</span>)</a></code></pre></div>
<p><img src="" /><!-- --></p>
<p>CPO值表明模型拟合较好，参考线以上的点的优势意味着对简化模型的偏好。这与我们之前使用DIC和LPML标准的研究结果一致。</p>
</div>
</div>
<div id="bayesian-residual-plots" class="section level2">
<h2>3.4 Bayesian residual plots</h2>
<p>使用<code>bri.lmresid.plot</code>生成贝叶斯残差图:</p>
<div class="sourceCode" id="cb35"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb35-1" title="1"><span class="kw">par</span>(<span class="dt">mfrow=</span><span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">2</span>))</a>
<a class="sourceLine" id="cb35-2" title="2"><span class="kw">bri.lmresid.plot</span>(usair.inla3)</a>
<a class="sourceLine" id="cb35-3" title="3"><span class="kw">bri.lmresid.plot</span>(usair.inla3, usair<span class="op">$</span>negtemp, <span class="dt">xlab =</span> <span class="st">&quot;negtemp&quot;</span>, <span class="dt">smooth =</span><span class="ot">TRUE</span>)</a></code></pre></div>
<p><img src="" /><!-- --></p>
<p>贝叶斯残差一般都在零附近呈现随机模式。但我们发现观测数31似乎是一个离群值，其贝叶斯残差高达59.6927。</p>
<hr />
</div>
</div>
<div id="robust-linear-regression-with-t-distribution" class="section level1">
<h1>4. Robust Linear regression with t-distribution</h1>
<p>在INLA中，通过指定<code>family=&quot;T&quot;</code>来实现。</p>
<div class="sourceCode" id="cb36"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb36-1" title="1"><span class="co">## Usair data example</span></a>
<a class="sourceLine" id="cb36-2" title="2">usair.inla4 &lt;-<span class="st"> </span><span class="kw">inla</span>(usair.formula2, <span class="dt">family =</span> <span class="st">&quot;T&quot;</span>, <span class="dt">data =</span> usair, <span class="dt">control.compute =</span> <span class="kw">list</span>(<span class="dt">dic =</span> <span class="ot">TRUE</span>, <span class="dt">cpo =</span> <span class="ot">TRUE</span>))</a>
<a class="sourceLine" id="cb36-3" title="3"><span class="co"># summary(usair.inla4)</span></a>
<a class="sourceLine" id="cb36-4" title="4"></a>
<a class="sourceLine" id="cb36-5" title="5">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">round</span>(usair.inla4<span class="op">$</span>summary.fixed, <span class="dv">4</span>),<span class="dt">caption =</span> <span class="st">&#39;固定效应情况&#39;</span>, <span class="dt">align=</span><span class="st">&#39;c&#39;</span>)</a></code></pre></div>
<table>
<caption>固定效应情况</caption>
<thead>
<tr class="header">
<th></th>
<th align="center">mean</th>
<th align="center">sd</th>
<th align="center">0.025quant</th>
<th align="center">0.5quant</th>
<th align="center">0.975quant</th>
<th align="center">mode</th>
<th align="center">kld</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>(Intercept)</td>
<td align="center">118.4034</td>
<td align="center">26.0576</td>
<td align="center">67.4027</td>
<td align="center">118.2333</td>
<td align="center">170.3176</td>
<td align="center">117.9043</td>
<td align="center">0</td>
</tr>
<tr class="even">
<td>negtemp</td>
<td align="center">1.4297</td>
<td align="center">0.3447</td>
<td align="center">0.7729</td>
<td align="center">1.4216</td>
<td align="center">2.1329</td>
<td align="center">1.4051</td>
<td align="center">0</td>
</tr>
<tr class="odd">
<td>manuf</td>
<td align="center">0.0267</td>
<td align="center">0.0036</td>
<td align="center">0.0192</td>
<td align="center">0.0268</td>
<td align="center">0.0335</td>
<td align="center">0.0270</td>
<td align="center">0</td>
</tr>
<tr class="even">
<td>wind</td>
<td align="center">-4.0148</td>
<td align="center">1.5948</td>
<td align="center">-7.1428</td>
<td align="center">-4.0229</td>
<td align="center">-0.8415</td>
<td align="center">-4.0363</td>
<td align="center">0</td>
</tr>
<tr class="odd">
<td>precip</td>
<td align="center">0.4257</td>
<td align="center">0.1909</td>
<td align="center">0.0682</td>
<td align="center">0.4186</td>
<td align="center">0.8223</td>
<td align="center">0.4038</td>
<td align="center">0</td>
</tr>
</tbody>
</table>
<div class="sourceCode" id="cb37"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb37-1" title="1">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">round</span>(usair.inla4<span class="op">$</span>summary.hyperpar, <span class="dv">4</span>), <span class="dt">caption =</span> <span class="st">&#39;</span></a>
<a class="sourceLine" id="cb37-2" title="2"><span class="st">    超参数情况&#39;</span>, <span class="dt">align=</span><span class="st">&#39;c&#39;</span>)</a></code></pre></div>
<table>
<caption>超参数情况</caption>
<thead>
<tr class="header">
<th></th>
<th align="center">mean</th>
<th align="center">sd</th>
<th align="center">0.025quant</th>
<th align="center">0.5quant</th>
<th align="center">0.975quant</th>
<th align="center">mode</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>precision for the student-t observations</td>
<td align="center">0.0053</td>
<td align="center">0.0016</td>
<td align="center">0.003</td>
<td align="center">0.0051</td>
<td align="center">0.0092</td>
<td align="center">0.0046</td>
</tr>
<tr class="even">
<td>degrees of freedom for student-t</td>
<td align="center">10.9806</td>
<td align="center">8.7815</td>
<td align="center">3.567</td>
<td align="center">8.3861</td>
<td align="center">34.1006</td>
<td align="center">5.6094</td>
</tr>
</tbody>
</table>
</div>
<div id="analysis-of-variance" class="section level1">
<h1>5. Analysis of Variance</h1>
<p>下面的例子，我们将只关注<strong>固定效应</strong>的方差分析模型。关于随机效应模型，将在第5章详细介绍。</p>
<p>来自一项研究可待因和针灸对男性患者术后牙痛的影响的实验(Kutner等，2004年)。</p>
<p>该研究采用随机区组设计，在一个因子结构中出现两个治疗因素。两种治疗因素都有两个层次。</p>
<p>反应变量缓解是疼痛缓解评分(分数越高，患者缓解程度越好)。根据对疼痛耐受性的评估，32名受试者被分配到8个区块，每个区块4名受试者。</p>
<div class="sourceCode" id="cb38"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb38-1" title="1"><span class="co"># PainRelief data</span></a>
<a class="sourceLine" id="cb38-2" title="2"><span class="kw">data</span>(painrelief, <span class="dt">package =</span> <span class="st">&quot;brinla&quot;</span>)</a>
<a class="sourceLine" id="cb38-3" title="3"></a>
<a class="sourceLine" id="cb38-4" title="4">painrelief<span class="op">$</span>PainLevel &lt;-<span class="st"> </span><span class="kw">as.factor</span>(painrelief<span class="op">$</span>PainLevel)</a>
<a class="sourceLine" id="cb38-5" title="5">painrelief<span class="op">$</span>Codeine &lt;-<span class="st"> </span><span class="kw">as.factor</span>(painrelief<span class="op">$</span>Codeine)</a>
<a class="sourceLine" id="cb38-6" title="6">painrelief<span class="op">$</span>Acupuncture &lt;-<span class="st"> </span><span class="kw">as.factor</span>(painrelief<span class="op">$</span>Acupuncture)</a>
<a class="sourceLine" id="cb38-7" title="7"></a>
<a class="sourceLine" id="cb38-8" title="8">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">head</span>(painrelief), <span class="dt">caption =</span> <span class="st">&#39;</span></a>
<a class="sourceLine" id="cb38-9" title="9"><span class="st">    painrelief数据&#39;</span>, <span class="dt">align=</span><span class="st">&#39;c&#39;</span>)</a></code></pre></div>
<table>
<caption>painrelief数据</caption>
<thead>
<tr class="header">
<th align="center">PainLevel</th>
<th align="center">Codeine</th>
<th align="center">Acupuncture</th>
<th align="center">Relief</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td align="center">1</td>
<td align="center">1</td>
<td align="center">1</td>
<td align="center">0.0</td>
</tr>
<tr class="even">
<td align="center">1</td>
<td align="center">2</td>
<td align="center">1</td>
<td align="center">0.5</td>
</tr>
<tr class="odd">
<td align="center">1</td>
<td align="center">1</td>
<td align="center">2</td>
<td align="center">0.6</td>
</tr>
<tr class="even">
<td align="center">1</td>
<td align="center">2</td>
<td align="center">2</td>
<td align="center">1.2</td>
</tr>
<tr class="odd">
<td align="center">2</td>
<td align="center">1</td>
<td align="center">1</td>
<td align="center">0.3</td>
</tr>
<tr class="even">
<td align="center">2</td>
<td align="center">2</td>
<td align="center">1</td>
<td align="center">0.6</td>
</tr>
</tbody>
</table>
<div class="sourceCode" id="cb39"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb39-1" title="1">painrelief.inla &lt;-<span class="st"> </span><span class="kw">inla</span>(Relief <span class="op">~</span><span class="st"> </span>PainLevel <span class="op">+</span><span class="st"> </span>Codeine<span class="op">*</span>Acupuncture, <span class="dt">data =</span> painrelief)</a>
<a class="sourceLine" id="cb39-2" title="2"><span class="co">#summary(painrelief.inla)</span></a>
<a class="sourceLine" id="cb39-3" title="3">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">round</span>(painrelief.inla<span class="op">$</span>summary.fixed, <span class="dv">4</span>))</a></code></pre></div>
<table>
<thead>
<tr class="header">
<th></th>
<th align="right">mean</th>
<th align="right">sd</th>
<th align="right">0.025quant</th>
<th align="right">0.5quant</th>
<th align="right">0.975quant</th>
<th align="right">mode</th>
<th align="right">kld</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>(Intercept)</td>
<td align="right">0.0188</td>
<td align="right">0.0702</td>
<td align="right">-0.1202</td>
<td align="right">0.0188</td>
<td align="right">0.1577</td>
<td align="right">0.0188</td>
<td align="right">1e-04</td>
</tr>
<tr class="even">
<td>PainLevel2</td>
<td align="right">0.1500</td>
<td align="right">0.0846</td>
<td align="right">-0.0176</td>
<td align="right">0.1500</td>
<td align="right">0.3175</td>
<td align="right">0.1500</td>
<td align="right">1e-04</td>
</tr>
<tr class="odd">
<td>PainLevel3</td>
<td align="right">0.3250</td>
<td align="right">0.0846</td>
<td align="right">0.1574</td>
<td align="right">0.3250</td>
<td align="right">0.4925</td>
<td align="right">0.3250</td>
<td align="right">1e-04</td>
</tr>
<tr class="even">
<td>PainLevel4</td>
<td align="right">0.3000</td>
<td align="right">0.0846</td>
<td align="right">0.1324</td>
<td align="right">0.3000</td>
<td align="right">0.4675</td>
<td align="right">0.3000</td>
<td align="right">1e-04</td>
</tr>
<tr class="odd">
<td>PainLevel5</td>
<td align="right">0.6750</td>
<td align="right">0.0846</td>
<td align="right">0.5074</td>
<td align="right">0.6750</td>
<td align="right">0.8425</td>
<td align="right">0.6750</td>
<td align="right">1e-04</td>
</tr>
<tr class="even">
<td>PainLevel6</td>
<td align="right">0.9750</td>
<td align="right">0.0846</td>
<td align="right">0.8074</td>
<td align="right">0.9750</td>
<td align="right">1.1425</td>
<td align="right">0.9750</td>
<td align="right">1e-04</td>
</tr>
<tr class="odd">
<td>PainLevel7</td>
<td align="right">1.0750</td>
<td align="right">0.0846</td>
<td align="right">0.9074</td>
<td align="right">1.0750</td>
<td align="right">1.2425</td>
<td align="right">1.0750</td>
<td align="right">1e-04</td>
</tr>
<tr class="even">
<td>PainLevel8</td>
<td align="right">1.1500</td>
<td align="right">0.0846</td>
<td align="right">0.9824</td>
<td align="right">1.1500</td>
<td align="right">1.3175</td>
<td align="right">1.1500</td>
<td align="right">1e-04</td>
</tr>
<tr class="odd">
<td>Codeine2</td>
<td align="right">0.4625</td>
<td align="right">0.0598</td>
<td align="right">0.3440</td>
<td align="right">0.4625</td>
<td align="right">0.5810</td>
<td align="right">0.4625</td>
<td align="right">1e-04</td>
</tr>
<tr class="even">
<td>Acupuncture2</td>
<td align="right">0.5750</td>
<td align="right">0.0598</td>
<td align="right">0.4565</td>
<td align="right">0.5750</td>
<td align="right">0.6935</td>
<td align="right">0.5750</td>
<td align="right">1e-04</td>
</tr>
<tr class="odd">
<td>Codeine2:Acupuncture2</td>
<td align="right">0.1500</td>
<td align="right">0.0846</td>
<td align="right">-0.0176</td>
<td align="right">0.1500</td>
<td align="right">0.3175</td>
<td align="right">0.1500</td>
<td align="right">1e-04</td>
</tr>
</tbody>
</table>
<p>治疗因素可待因和针灸的主要作用在贝叶斯意义上都是高度显著的。但在95%可信水平上，两者之间的交互作用不显著，说明两者之间不存在交互作用。</p>
<div class="sourceCode" id="cb40"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb40-1" title="1">est1 &lt;-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">x =</span><span class="kw">c</span>(<span class="st">&quot;2&quot;</span>, <span class="st">&quot;3&quot;</span>, <span class="st">&quot;4&quot;</span>, <span class="st">&quot;5&quot;</span>, <span class="st">&quot;6&quot;</span>, <span class="st">&quot;7&quot;</span>, <span class="st">&quot;8&quot;</span>),</a>
<a class="sourceLine" id="cb40-2" title="2">                   <span class="dt">Estimate =</span> painrelief.inla<span class="op">$</span>summary.fixed[<span class="kw">c</span>(<span class="dv">2</span><span class="op">:</span><span class="dv">8</span>),<span class="dv">1</span>],</a>
<a class="sourceLine" id="cb40-3" title="3">                   <span class="dt">L =</span> painrelief.inla<span class="op">$</span>summary.fixed[<span class="kw">c</span>(<span class="dv">2</span><span class="op">:</span><span class="dv">8</span>),<span class="dv">3</span>],</a>
<a class="sourceLine" id="cb40-4" title="4">                   <span class="dt">U =</span> painrelief.inla<span class="op">$</span>summary.fixed[<span class="kw">c</span>(<span class="dv">2</span><span class="op">:</span><span class="dv">8</span>),<span class="dv">5</span>])</a></code></pre></div>
<p>为了更好地理解影响的重要性，我们生成了不同疼痛水平的后验均值估计和95%可信水平的图。水平线是疼痛级别1的参考线，它被设置为0。疼痛程度越高，受试者的疼痛缓解评分就越高(疼痛程度4除外)。显然，疼痛程度是模型中需要考虑的一个显著的混杂因素。</p>
<div class="sourceCode" id="cb41"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb41-1" title="1">p1 &lt;-<span class="st"> </span><span class="kw">ggplot</span>(est1, <span class="kw">aes</span>(<span class="dt">x =</span> x, <span class="dt">y =</span> Estimate)) <span class="op">+</span><span class="st"> </span><span class="kw">geom_point</span>(<span class="dt">size =</span> <span class="dv">5</span>) <span class="op">+</span><span class="st"> </span><span class="kw">geom_errorbar</span>(<span class="kw">aes</span>(<span class="dt">ymax =</span> U, <span class="dt">ymin =</span> L),<span class="dt">size=</span><span class="dv">1</span>) <span class="op">+</span><span class="st"> </span><span class="kw">geom_hline</span>(<span class="dt">yintercept=</span><span class="dv">0</span>, <span class="dt">size=</span><span class="fl">1.5</span>, <span class="dt">col=</span><span class="st">&quot;black&quot;</span>) <span class="op">+</span><span class="st"> </span><span class="kw">xlab</span>(<span class="st">&quot;Pain Level&quot;</span>)</a>
<a class="sourceLine" id="cb41-2" title="2">p1</a></code></pre></div>
<p><img src="" /><!-- --></p>
</div>
<div id="ridge-regression" class="section level1">
<h1>6. Ridge regression</h1>
<p><strong>法国经济进口活动有关数据</strong>：因变量为<code>进口(import)</code>、<code>国内生产(DOPROD)</code>、<code>股票形式(stock)</code>和<code>国内消费(CONSUM)</code>。</p>
<ul>
<li>样本相关性</li>
</ul>
<div class="sourceCode" id="cb42"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb42-1" title="1"><span class="kw">data</span>(frencheconomy, <span class="dt">package =</span> <span class="st">&quot;brinla&quot;</span>)</a>
<a class="sourceLine" id="cb42-2" title="2"><span class="kw">head</span>(frencheconomy)</a></code></pre></div>
<pre><code>##   YEAR IMPORT DOPROD STOCK CONSUM
## 1   49   15.9  149.3   4.2  108.1
## 2   50   16.4  161.2   4.1  114.8
## 3   51   19.0  171.5   3.1  123.2
## 4   52   19.1  175.5   3.1  126.9
## 5   53   18.8  180.8   1.1  132.1
## 6   54   20.4  190.7   2.2  137.7</code></pre>
<div class="sourceCode" id="cb44"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb44-1" title="1">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">round</span>(<span class="kw">cor</span>(frencheconomy[,<span class="op">-</span><span class="dv">1</span>]),<span class="dv">4</span>), <span class="dt">caption =</span> <span class="st">&#39;</span></a>
<a class="sourceLine" id="cb44-2" title="2"><span class="st">    frencheconomy数据&#39;</span>, <span class="dt">align=</span><span class="st">&#39;c&#39;</span>)</a></code></pre></div>
<table>
<caption>frencheconomy数据</caption>
<thead>
<tr class="header">
<th></th>
<th align="center">IMPORT</th>
<th align="center">DOPROD</th>
<th align="center">STOCK</th>
<th align="center">CONSUM</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>IMPORT</td>
<td align="center">1.0000</td>
<td align="center">0.9842</td>
<td align="center">0.2659</td>
<td align="center">0.9848</td>
</tr>
<tr class="even">
<td>DOPROD</td>
<td align="center">0.9842</td>
<td align="center">1.0000</td>
<td align="center">0.2154</td>
<td align="center">0.9989</td>
</tr>
<tr class="odd">
<td>STOCK</td>
<td align="center">0.2659</td>
<td align="center">0.2154</td>
<td align="center">1.0000</td>
<td align="center">0.2137</td>
</tr>
<tr class="even">
<td>CONSUM</td>
<td align="center">0.9848</td>
<td align="center">0.9989</td>
<td align="center">0.2137</td>
<td align="center">1.0000</td>
</tr>
</tbody>
</table>
<ul>
<li>数据标准化</li>
</ul>
<p>贝叶斯岭回归假设所有预测因子的系数元素(<span class="math inline">\(\beta_1,\dots,\beta_p\)</span>)都是从一个标准正态密度中提取的。</p>
<div class="sourceCode" id="cb45"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb45-1" title="1">fe.scaled &lt;-<span class="st"> </span><span class="kw">cbind</span>(frencheconomy[, <span class="dv">1</span><span class="op">:</span><span class="dv">2</span>], <span class="kw">scale</span>(frencheconomy[, <span class="kw">c</span>(<span class="op">-</span><span class="dv">1</span>,<span class="op">-</span><span class="dv">2</span>)]))</a></code></pre></div>
<div class="sourceCode" id="cb46"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb46-1" title="1"><span class="co"># set priors</span></a>
<a class="sourceLine" id="cb46-2" title="2">n &lt;-<span class="st"> </span><span class="kw">nrow</span>(frencheconomy)</a>
<a class="sourceLine" id="cb46-3" title="3"></a>
<a class="sourceLine" id="cb46-4" title="4">fe.scaled<span class="op">$</span>beta1 &lt;-<span class="st"> </span><span class="kw">rep</span>(<span class="dv">1</span>,n)</a>
<a class="sourceLine" id="cb46-5" title="5">fe.scaled<span class="op">$</span>beta2 &lt;-<span class="st"> </span><span class="kw">rep</span>(<span class="dv">2</span>,n)</a>
<a class="sourceLine" id="cb46-6" title="6">fe.scaled<span class="op">$</span>beta3 &lt;-<span class="st"> </span><span class="kw">rep</span>(<span class="dv">3</span>,n)</a></code></pre></div>
<p>对于岭回归，参数的先验有一个共同的未知方差，要实现这一点，我们必须使用副本，并改变数据集</p>
<div class="sourceCode" id="cb47"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb47-1" title="1"><span class="co"># this is the prior for the precision of beta</span></a>
<a class="sourceLine" id="cb47-2" title="2">param.beta =<span class="st"> </span><span class="kw">list</span>(<span class="dt">prec =</span> <span class="kw">list</span>(<span class="dt">param =</span> <span class="kw">c</span>(<span class="fl">1.0e-3</span>, <span class="fl">1.0e-3</span>)))</a>
<a class="sourceLine" id="cb47-3" title="3"></a>
<a class="sourceLine" id="cb47-4" title="4">formula.ridge  =<span class="st"> </span>IMPORT <span class="op">~</span><span class="st"> </span><span class="kw">f</span>(beta1, DOPROD, <span class="dt">model=</span><span class="st">&quot;iid&quot;</span>, <span class="dt">values =</span> <span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">2</span>,<span class="dv">3</span>), <span class="dt">hyper =</span> param.beta) <span class="op">+</span><span class="st"> </span><span class="kw">f</span>(beta2, STOCK, <span class="dt">copy=</span><span class="st">&quot;beta1&quot;</span>, <span class="dt">fixed=</span>T) <span class="op">+</span><span class="st"> </span><span class="kw">f</span>(beta3, CONSUM, <span class="dt">copy=</span><span class="st">&quot;beta1&quot;</span>, <span class="dt">fixed=</span>T)</a>
<a class="sourceLine" id="cb47-5" title="5">frencheconomy.ridge &lt;-<span class="st"> </span><span class="kw">inla</span>(formula.ridge, <span class="dt">data=</span>fe.scaled)</a>
<a class="sourceLine" id="cb47-6" title="6">ridge.est &lt;-<span class="st"> </span><span class="kw">rbind</span>(frencheconomy.ridge<span class="op">$</span>summary.fixed, frencheconomy.ridge<span class="op">$</span>summary.random<span class="op">$</span>beta1[,<span class="op">-</span><span class="dv">1</span>]) </a>
<a class="sourceLine" id="cb47-7" title="7"></a>
<a class="sourceLine" id="cb47-8" title="8">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">round</span>(ridge.est,<span class="dv">4</span>))</a></code></pre></div>
<table>
<thead>
<tr class="header">
<th></th>
<th align="right">mean</th>
<th align="right">sd</th>
<th align="right">0.025quant</th>
<th align="right">0.5quant</th>
<th align="right">0.975quant</th>
<th align="right">mode</th>
<th align="right">kld</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>(Intercept)</td>
<td align="right">30.0778</td>
<td align="right">0.5198</td>
<td align="right">29.0449</td>
<td align="right">30.0778</td>
<td align="right">31.109</td>
<td align="right">30.0778</td>
<td align="right">0e+00</td>
</tr>
<tr class="even">
<td>1</td>
<td align="right">5.1131</td>
<td align="right">5.4237</td>
<td align="right">-6.7892</td>
<td align="right">5.3431</td>
<td align="right">15.668</td>
<td align="right">5.6205</td>
<td align="right">3e-04</td>
</tr>
<tr class="odd">
<td>2</td>
<td align="right">0.7205</td>
<td align="right">0.5451</td>
<td align="right">-0.3618</td>
<td align="right">0.7202</td>
<td align="right">1.802</td>
<td align="right">0.7199</td>
<td align="right">0e+00</td>
</tr>
<tr class="even">
<td>3</td>
<td align="right">6.9769</td>
<td align="right">5.4265</td>
<td align="right">-3.5327</td>
<td align="right">6.7312</td>
<td align="right">18.923</td>
<td align="right">6.4159</td>
<td align="right">3e-04</td>
</tr>
</tbody>
</table>
<ul>
<li>Comparing with Standard Bayesian Linear regression</li>
</ul>
<div class="sourceCode" id="cb48"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb48-1" title="1">formula &lt;-<span class="st"> </span>IMPORT <span class="op">~</span><span class="st">  </span>DOPROD <span class="op">+</span><span class="st"> </span>STOCK <span class="op">+</span><span class="st"> </span>CONSUM</a>
<a class="sourceLine" id="cb48-2" title="2">frencheconomy.inla &lt;-<span class="st"> </span><span class="kw">inla</span>(formula, <span class="dt">data =</span> fe.scaled, <span class="dt">control.fixed =</span> <span class="kw">list</span>(<span class="dt">prec =</span> <span class="fl">1.0e-3</span>), <span class="dt">control.family =</span> <span class="kw">list</span>(<span class="dt">hyper =</span> param.beta))</a>
<a class="sourceLine" id="cb48-3" title="3">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">round</span>(frencheconomy.inla<span class="op">$</span>summary.fixed, <span class="dv">4</span>))</a></code></pre></div>
<table>
<thead>
<tr class="header">
<th></th>
<th align="right">mean</th>
<th align="right">sd</th>
<th align="right">0.025quant</th>
<th align="right">0.5quant</th>
<th align="right">0.975quant</th>
<th align="right">mode</th>
<th align="right">kld</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>(Intercept)</td>
<td align="right">30.0778</td>
<td align="right">0.5729</td>
<td align="right">28.9378</td>
<td align="right">30.0778</td>
<td align="right">31.216</td>
<td align="right">30.0778</td>
<td align="right">0</td>
</tr>
<tr class="even">
<td>DOPROD</td>
<td align="right">3.0052</td>
<td align="right">10.9358</td>
<td align="right">-18.5325</td>
<td align="right">2.9662</td>
<td align="right">24.743</td>
<td align="right">2.8961</td>
<td align="right">0</td>
</tr>
<tr class="odd">
<td>STOCK</td>
<td align="right">0.7197</td>
<td align="right">0.6037</td>
<td align="right">-0.4819</td>
<td align="right">0.7198</td>
<td align="right">1.919</td>
<td align="right">0.7199</td>
<td align="right">0</td>
</tr>
<tr class="even">
<td>CONSUM</td>
<td align="right">9.1322</td>
<td align="right">10.9315</td>
<td align="right">-12.6218</td>
<td align="right">9.1707</td>
<td align="right">30.635</td>
<td align="right">9.2416</td>
<td align="right">0</td>
</tr>
</tbody>
</table>
<ul>
<li>Comparing with Ridge regression frequentist approach</li>
</ul>
<div class="sourceCode" id="cb49"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb49-1" title="1">reg2 &lt;-<span class="st"> </span><span class="kw">lm.ridge</span>(IMPORT <span class="op">~</span><span class="st">  </span>DOPROD <span class="op">+</span><span class="st"> </span>STOCK <span class="op">+</span><span class="st"> </span>CONSUM, <span class="dt">data =</span> fe.scaled, <span class="dt">lambda =</span> <span class="kw">seq</span>(<span class="dv">0</span>, <span class="dv">1</span>, <span class="dt">length=</span><span class="dv">100</span>))</a>
<a class="sourceLine" id="cb49-2" title="2"></a>
<a class="sourceLine" id="cb49-3" title="3">reg2.final &lt;-<span class="st"> </span><span class="kw">lm.ridge</span>(IMPORT <span class="op">~</span><span class="st">  </span>DOPROD <span class="op">+</span><span class="st"> </span>STOCK <span class="op">+</span><span class="st"> </span>CONSUM, <span class="dt">data =</span> fe.scaled, <span class="dt">lambda =</span> reg2<span class="op">$</span>kHKB)</a>
<a class="sourceLine" id="cb49-4" title="4">reg2.final</a></code></pre></div>
<pre><code>##          DOPROD   STOCK  CONSUM 
## 30.0778  4.9560  0.7176  7.1669</code></pre>
</div>
<div id="linear-regression-with-autoregressive-errors" class="section level1">
<h1>7. Linear regression with autoregressive errors</h1>
<p><strong>时间序列数据</strong>：新西兰的失业数据包括青年(15-19岁)和成人(19岁以上)的季度失业率。</p>
<p>自2008年6月起，新西兰政府废除了《最低工资法》。在此，我们想研究在该法案废除之前和之后，成年人和年轻人失业率之间的关系。</p>
<div class="sourceCode" id="cb51"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb51-1" title="1"><span class="co"># Read the data</span></a>
<a class="sourceLine" id="cb51-2" title="2"><span class="kw">data</span>(nzunemploy, <span class="dt">package =</span> <span class="st">&quot;brinla&quot;</span>)</a>
<a class="sourceLine" id="cb51-3" title="3">nzunemploy<span class="op">$</span>time &lt;-<span class="st"> </span><span class="dv">1</span><span class="op">:</span><span class="kw">nrow</span>(nzunemploy)</a></code></pre></div>
<p>绘制成人和青年的时间序列数据</p>
<div class="sourceCode" id="cb52"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb52-1" title="1"><span class="co">#library(tidyr)</span></a>
<a class="sourceLine" id="cb52-2" title="2"><span class="kw">qplot</span>(time, value, <span class="dt">data =</span> <span class="kw">gather</span>(nzunemploy[,<span class="kw">c</span>(<span class="dv">2</span>,<span class="dv">3</span>,<span class="dv">5</span>)], variable, value, <span class="op">-</span>time), <span class="dt">geom =</span> <span class="st">&quot;line&quot;</span>) <span class="op">+</span><span class="st"> </span><span class="kw">geom_vline</span>(<span class="dt">xintercept =</span> <span class="dv">90</span>) <span class="op">+</span><span class="st"> </span><span class="kw">facet_grid</span>(variable <span class="op">~</span><span class="st"> </span>., <span class="dt">scale =</span> <span class="st">&quot;free&quot;</span>) <span class="op">+</span><span class="st"> </span><span class="kw">ylab</span>(<span class="st">&quot;Unemployment rate&quot;</span>) <span class="op">+</span><span class="st"> </span><span class="kw">theme_bw</span>()</a></code></pre></div>
<p><img src="" /><!-- --></p>
<p>为了使系数更容易理解，我们将成人失业率集中在时间序列的平均值上。</p>
<div class="sourceCode" id="cb53"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb53-1" title="1"><span class="co"># Centering predictor</span></a>
<a class="sourceLine" id="cb53-2" title="2">nzunemploy<span class="op">$</span>centeredadult =<span class="st"> </span><span class="kw">with</span>(nzunemploy, adult <span class="op">-</span><span class="st"> </span><span class="kw">mean</span>(adult))</a></code></pre></div>
<ul>
<li>拟合一个具有独立误差的标准线性回归</li>
</ul>
<div class="sourceCode" id="cb54"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb54-1" title="1">formula1 &lt;-<span class="st"> </span>youth <span class="op">~</span><span class="st"> </span>centeredadult<span class="op">*</span>policy</a>
<a class="sourceLine" id="cb54-2" title="2">nzunemploy.inla1 &lt;-<span class="st"> </span><span class="kw">inla</span>(formula1, <span class="dt">data=</span> nzunemploy)</a>
<a class="sourceLine" id="cb54-3" title="3"><span class="co"># summary(nzunemploy.inla1)</span></a>
<a class="sourceLine" id="cb54-4" title="4"><span class="kw">round</span>(nzunemploy.inla1<span class="op">$</span>summary.fixed, <span class="dv">4</span>)</a></code></pre></div>
<pre><code>##                             mean     sd 0.025quant 0.5quant 0.975quant   mode
## (Intercept)               16.282 0.1534     15.980   16.282     16.584 16.282
## centeredadult              1.533 0.0750      1.386    1.533      1.681  1.533
## policyEqual                9.442 0.5258      8.407    9.442     10.475  9.442
## centeredadult:policyEqual  2.853 0.4616      1.945    2.853      3.761  2.853
##                           kld
## (Intercept)                 0
## centeredadult               0
## policyEqual                 0
## centeredadult:policyEqual   0</code></pre>
<div class="sourceCode" id="cb56"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb56-1" title="1"><span class="co"># Plot the Bayesian residuals</span></a>
<a class="sourceLine" id="cb56-2" title="2"><span class="kw">par</span>(<span class="dt">mfrow=</span><span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">1</span>))</a>
<a class="sourceLine" id="cb56-3" title="3">nzunemploy.res1 &lt;-<span class="st"> </span><span class="kw">bri.lmresid.plot</span>(nzunemploy.inla1, <span class="dt">type=</span><span class="st">&quot;o&quot;</span>)</a></code></pre></div>
<p><img src="" /><!-- --></p>
<p>它们之间存在一定程度的自相关。</p>
<div class="sourceCode" id="cb57"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb57-1" title="1"><span class="co"># Plot the autocorrelation and partial autocorrelation</span></a>
<a class="sourceLine" id="cb57-2" title="2"><span class="kw">par</span>(<span class="dt">mfrow =</span> <span class="kw">c</span>(<span class="dv">2</span>,<span class="dv">1</span>))</a>
<a class="sourceLine" id="cb57-3" title="3"><span class="kw">acf</span>(nzunemploy.res1<span class="op">$</span>resid, <span class="dt">main =</span> <span class="st">&quot;&quot;</span>)</a>
<a class="sourceLine" id="cb57-4" title="4"><span class="kw">acf</span>(nzunemploy.res1<span class="op">$</span>resid, <span class="dt">type =</span> <span class="st">&quot;partial&quot;</span>, <span class="dt">main=</span><span class="st">&quot;&quot;</span>)</a></code></pre></div>
<p><img src="" /><!-- --></p>
<p>图上的虚线对应95%置信带。自相关函数的形式呈指数衰减，而偏自相关函数的形式在滞后1处有很高的峰值。这些表明AR(1)过程将适合于回归模型中的误差项。</p>
<ul>
<li>拟合一个有AR(1)误差的线性回归</li>
</ul>
<div class="sourceCode" id="cb58"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb58-1" title="1">formula2 &lt;-<span class="st"> </span>youth <span class="op">~</span><span class="st"> </span>centeredadult<span class="op">*</span>policy <span class="op">+</span><span class="st"> </span><span class="kw">f</span>(time, <span class="dt">model =</span> <span class="st">&quot;ar1&quot;</span>)</a>
<a class="sourceLine" id="cb58-2" title="2">nzunemploy.inla2 &lt;-<span class="st"> </span><span class="kw">inla</span>(formula2, <span class="dt">data =</span> nzunemploy, <span class="dt">control.family =</span> <span class="kw">list</span>(<span class="dt">hyper =</span> <span class="kw">list</span>(<span class="dt">prec =</span> <span class="kw">list</span>(<span class="dt">initial =</span> <span class="dv">15</span>, <span class="dt">fixed =</span> <span class="ot">TRUE</span>))))</a>
<a class="sourceLine" id="cb58-3" title="3"><span class="co"># summary(nzunemploy.inla2)</span></a>
<a class="sourceLine" id="cb58-4" title="4">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">round</span>(nzunemploy.inla2<span class="op">$</span>summary.fixed, <span class="dv">4</span>))</a></code></pre></div>
<table>
<thead>
<tr class="header">
<th></th>
<th align="right">mean</th>
<th align="right">sd</th>
<th align="right">0.025quant</th>
<th align="right">0.5quant</th>
<th align="right">0.975quant</th>
<th align="right">mode</th>
<th align="right">kld</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>(Intercept)</td>
<td align="right">16.344</td>
<td align="right">0.3029</td>
<td align="right">15.768</td>
<td align="right">16.334</td>
<td align="right">16.975</td>
<td align="right">16.322</td>
<td align="right">5e-04</td>
</tr>
<tr class="even">
<td>centeredadult</td>
<td align="right">1.520</td>
<td align="right">0.1354</td>
<td align="right">1.246</td>
<td align="right">1.521</td>
<td align="right">1.785</td>
<td align="right">1.524</td>
<td align="right">1e-04</td>
</tr>
<tr class="odd">
<td>policyEqual</td>
<td align="right">8.981</td>
<td align="right">0.9724</td>
<td align="right">6.876</td>
<td align="right">9.036</td>
<td align="right">10.746</td>
<td align="right">9.117</td>
<td align="right">8e-04</td>
</tr>
<tr class="even">
<td>centeredadult:policyEqual</td>
<td align="right">2.516</td>
<td align="right">0.6009</td>
<td align="right">1.302</td>
<td align="right">2.525</td>
<td align="right">3.672</td>
<td align="right">2.543</td>
<td align="right">1e-04</td>
</tr>
</tbody>
</table>
<div class="sourceCode" id="cb59"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb59-1" title="1">knitr<span class="op">::</span><span class="kw">kable</span>(<span class="kw">round</span>(nzunemploy.inla2<span class="op">$</span>summary.hyperpar, <span class="dv">4</span>))</a></code></pre></div>
<table>
<thead>
<tr class="header">
<th></th>
<th align="right">mean</th>
<th align="right">sd</th>
<th align="right">0.025quant</th>
<th align="right">0.5quant</th>
<th align="right">0.975quant</th>
<th align="right">mode</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>Precision for time</td>
<td align="right">0.4464</td>
<td align="right">0.0859</td>
<td align="right">0.293</td>
<td align="right">0.4420</td>
<td align="right">0.6282</td>
<td align="right">0.4353</td>
</tr>
<tr class="even">
<td>Rho for time</td>
<td align="right">0.5128</td>
<td align="right">0.0898</td>
<td align="right">0.329</td>
<td align="right">0.5151</td>
<td align="right">0.6799</td>
<td align="right">0.5163</td>
</tr>
</tbody>
</table>
<p>注意，回归模型的精度prec固定在<span class="math inline">\(\tau = \exp(15)\)</span></p>
<ul>
<li>与频率派方法进行比对</li>
</ul>
<div class="sourceCode" id="cb60"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb60-1" title="1"><span class="co">#library(nlme)</span></a>
<a class="sourceLine" id="cb60-2" title="2">nzunemploy.gls &lt;-<span class="st"> </span><span class="kw">gls</span>(youth <span class="op">~</span><span class="st"> </span>centeredadult<span class="op">*</span>policy, <span class="dt">correlation =</span> <span class="kw">corAR1</span>(<span class="dt">form=</span><span class="op">~</span><span class="dv">1</span>), <span class="dt">data =</span> nzunemploy)</a>
<a class="sourceLine" id="cb60-3" title="3"><span class="kw">summary</span>(nzunemploy.gls)</a></code></pre></div>
<pre><code>## Generalized least squares fit by REML
##   Model: youth ~ centeredadult * policy 
##   Data: nzunemploy 
##   AIC   BIC logLik
##   353 368.5 -170.5
## 
## Correlation Structure: AR(1)
##  Formula: ~1 
##  Parameter estimate(s):
##    Phi 
## 0.5012 
## 
## Coefficients:
##                            Value Std.Error t-value p-value
## (Intercept)               16.329    0.2733   59.74       0
## centeredadult              1.522    0.1274   11.94       0
## policyEqual                9.083    0.8614   10.54       0
## centeredadult:policyEqual  2.545    0.5772    4.41       0
## 
##  Correlation: 
##                           (Intr) cntrdd plcyEq
## centeredadult             -0.020              
## policyEqual               -0.318  0.007       
## centeredadult:policyEqual -0.067 -0.155  0.583
## 
## Standardized residuals:
##     Min      Q1     Med      Q3     Max 
## -2.8923 -0.5546 -0.0242  0.5545  2.2957 
## 
## Residual standard error: 1.505 
## Degrees of freedom: 102 total; 98 residual</code></pre>
<div class="sourceCode" id="cb62"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb62-1" title="1"><span class="co"># plot the fitted lines</span></a>
<a class="sourceLine" id="cb62-2" title="2"><span class="kw">ggplot</span>(nzunemploy, <span class="kw">aes</span>(centeredadult, youth)) <span class="op">+</span><span class="st"> </span><span class="kw">geom_point</span>(<span class="kw">aes</span>(<span class="dt">shape =</span> <span class="kw">factor</span>(policy)), <span class="dt">size =</span> <span class="dv">3</span>) <span class="op">+</span><span class="st"> </span><span class="kw">geom_abline</span>(<span class="dt">intercept =</span> nzunemploy.inla2<span class="op">$</span>summary.fixed<span class="op">$</span>mean[<span class="dv">1</span>], <span class="dt">slope =</span> nzunemploy.inla2<span class="op">$</span>summary.fixed<span class="op">$</span>mean[<span class="dv">2</span>]) <span class="op">+</span><span class="st"> </span><span class="kw">geom_abline</span>(<span class="dt">intercept =</span> nzunemploy.inla2<span class="op">$</span>summary.fixed<span class="op">$</span>mean[<span class="dv">1</span>] <span class="op">+</span><span class="st"> </span>nzunemploy.inla2<span class="op">$</span>summary.fixed<span class="op">$</span>mean[<span class="dv">3</span>], <span class="dt">slope =</span> nzunemploy.inla2<span class="op">$</span>summary.fixed<span class="op">$</span>mean[<span class="dv">2</span>] <span class="op">+</span><span class="st"> </span>nzunemploy.inla2<span class="op">$</span>summary.fixed<span class="op">$</span>mean[<span class="dv">4</span>]) <span class="op">+</span><span class="st"> </span><span class="kw">theme_bw</span>()</a></code></pre></div>
<p><img src="" /><!-- --></p>
</div>
</section>



<!-- code folding -->


<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
  (function () {
    var script = document.createElement("script");
    script.type = "text/javascript";
    script.src  = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
    document.getElementsByTagName("head")[0].appendChild(script);
  })();
</script>

</body>
</html>
